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SUMMARY 


The objective of this study is to provide the numerical sim- 
ulation of the transsonic flows of idealized fluids and of incom- 
pressible viscous fluids, by the non linear least squares methods 
of R. GELOWINSKI and 0. PIRONNEAU. The complexity of the geometries 
studied in industrial aerodynamics explains the preference given to 
the finite elements for the approximation of the equations. 

Chapters 1, 2, 3, ^ describe the non linear equations, the 
boundary conditions and the various constraints controlling the 
two types of flow. The standard iterative methods for solving a 
quasi elliptical non linear equation with partial derivatives (E.D. 
P. ) are briefly reviewed in Chapter 5 with emphasis placed on two 
examples ; the fixed point method applied to the Gelder functional 
in the case of compressible subsonic flows and the Newton method 
used in the technique of decomposition of the lifting potential. 

Chapter 6 presents the new abstract least squares method. It 
consists of substituting the non linear equation by a problem of 
minimization in a H”"-*- type Sobolev functional space, which is itself 
equivalent to an optimal control problem and solved by a conjugate 
gradient algorithm with metric H 1 . The application of this method- 
ology to transsonic equations is presented in Chapter 7» We show 
how to include within the optimal control formulation two con- 
straints of aerodynamics; the condition of entropy, on the one hand, 
treated either by penalization or by artificial viscosity , and the 
Joukowski condition, on the other hand, taken into account by a fix- 
ed point method on circulation. 

The Navier-Stokes equations are reduced to a problem of minimi- 
zation in in the same manner in Chapter 8. Accordingly, we 

show that the state systems of the mixed optimal control problem 
are generalized Stolces problems in steady and unsteady cases, after 
quantification in time with the use of implicit Crank-Nicholson 
(for example) type schemes. To solve them, a mixed formulation 
proposed by GLOVINSKI-PIRONNEAU and based on certain decomposition 
properties of the biharmonic operator, is used. The Stokes algo- 
rithm is substitued by a sequence of Dirichlet problems coupled 
with an integral equation (e) conditioned on the pressure trace . 
defined on the boundary of the domain occupied by the fluid. 

Chapters 9 and 10 are devoted to the approximation of a trans- 
sonic and Navier-Stokes optimal control formulation by Pj- Lagrange 
conform finite elements, with degree lc=l or 2. The numerical imple- 
mentation of the conjugate gradient algorithms is developed and 
presented in the form of flow charts. Tho numerical implementa- 
tion of the Stokes algorithm ( Kj ) is described and the choice of a 
direct (Choloski) or iterative (preconditioned conjugate gradient) 
method for solving it is discussed. 

The large amounts of computations, due to complex tridimen- 
sional configurations (nacelle, vehicle, air-inlet, airplane) , 


stored in the main core of the computer* require an incomplete 

Choloski factorization of the discrete Dirichlet matrices shown 
on the inside of the control loop. The use of auxiliary operators 
LL C the solution of an optimal control problem is presented in 

Chapter 11 through comparisons of research results of J.A, MEIJ- 
ERFINK-M.A. VAN DER VORST and 0. AXELSSON. 

The numerical experiments are described in Chapter 12. The 
transsonic calculations obtained from the finite elements-optimal 
control codes are compared with those obtained from the finite 
differences codes of A. JAMESON on a NACA C012 airfoil and a Korn 
airfoil. 

More complex transsonic configurations of industrial aerody- 
namics such as multi-bodies or air inlets a* o analyzed. 

The feasibility of optimal control conjugate gradient algo- 
rithms is verified on bi and tridimensional Navier-Stokes calcu- 
lations, requiring considerable data processing resources (memory 
and CPU), Separated flows around/in an air inlet and around an 
swept-back wing with high incidence, are simulated numerically 
by following at various time cycles the evolution of the field of 
velocities, the field of pressures, the streamlines and the vorti- 
city. 


Finally, the last paragraph of Chapter 12 is devoted to the 
data processing efficiency of the auxiliary operators . It shows, 
through examples taken from the two flow families, how it is pos- 
sible, by using preconditioned optimal control algorithms, to cal- 
culate entirely in the main core of the computer, with small per- 
centages of Dirichlet matrices 7 , c without reducing the 

A d/I00 (5 - d<20) 

convergence velocity of the algorithm. 
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SOLUTION OF A FEW NON LINEAR PROBLEMS IN AERODYNAMICS BY THE FINITE 
ELEMENTS AND FUNCTIONAL LEAST SQUARES METHODS 

Jacques Periaux 

Pierre and Marie Curie University 


O. INTRODUCTION 

0.1. APPLICATIONS OF NON LINEAR AERODYNAMICS TO THE AERONAUTICS 
INDUSTRY 


The calculation of pressures in aeronautics plays an essential 
role in the optimization of aerodynamic shapes. The appearance of 
more and more powerful computers , over the past decade , both with 
respect to calculation speed and to memory capacity, has made it 
possible for the aviator to simulate numerically flows which ap- 
proximate more and more the flight conditions. To accomplish this 
it was necessary to define theoretically and numerically two fam- 
ilies of non linear equations : irrotational compressible idealized 
fluids, on the one hand, in order to study the transonic domain of 
the airplane, and incompressible viscuous fluids, on the other hand, 
modeled by the Navier-Stokes equations to provide a robust tool re- 
quired for the study of separated laminar flows in a first phase, 
then of turbulent flows in a second phase. 

The domaines occupied by the fluid are bi and tridimensional. 
They belong either to external aerodynamics when relating to air- 
foils (p) or to wings IVJ, or to internal aerodynamics when rela- 
ting to pipes (T), cavities (CA) or conduits (C). Finally, air in- 
lets belong to a third category ; mixed aerodynamics. The common 
denominator of these domaines is the complexity of the boundaries 
(one region surrounding a multibody, or one 3-U air inlrt composed 
of extremely complicated geometries, making it difficult to reduce 
it by conform conversion to a standard rectangular (or cubical) do- 
main !)• Furthermore, the final selection of the physical space as 
calculation domain was subjected to a numerical method by taking 
into account the boundary conditions with fine accuracy : THE FINITE 
ELEMENTS . 


0.2. Difficulties with respect to industrial configurations 

The numerical analysis of flows around industrial obstacles 
points up 3 types of difficulties t 

1 - geometrical difficulties! the configurations studied are 
extremely complex and require a delicate collection of data (descr- 
iption of a 2-D multi-body, or a wing + fuselage + air inlet + em- 
pennage type airplane configuration). 

2 - theoretical difficulties : the equations to be solved are 
non linear and their solutions may be composed of discontinuities. 
Furthermore, the following constraints must be satisfied simultan- 
eously x 




•constraint of aerodynamic reaction or Joukowski condition 
for perfect fluids , 

•constraint of physical shock or condition of entropy for 
transonic perfect fluids, 

•constraint of incompressiblity for viscuous fluids* 

3 - numerical difficulties t the volume of tridimensional calcu- 
lations (several thousands of unknowns) make it necessary to use al- 
gorithms which are both rapid for convergence and robust for stabil- 
ity. 


Figure 1 summarizes the situation in industry and describes the 
solution selected. 


1. - FEASIBLE MODELING OF AN INCOMPRESSIBLE IDEALIZED FLOW 
1,1, 2-D Non Lifting Case 

If ft and f designate respectively the domain and boundary 
of the region occupied by the fluid, as the latter is incompres- 
sible and irrotational , it obeys the following equations and boun- 
dary conditions (l) 


7*u ■ 0 ; u continuous 

7au - 0 (l) 

- g <r>, r-r^uTp 


where, in (l), u designates the fluid velocity and g the normal 
component of velocity on r ; on boundary sufficiently removed from 
the obstacle to ensure that the latter does not perturb the flow, 
g ■ u .JJ 1 where is the external standard of the domain, whereas 
on "wall of the obstacle (P) g=0 and ■* — ■ n means then that the 

r U*n 

fluid slides over the wall. 
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TIONED CONJUGATE GRADIENT 



The irrotational condition is expressed in a standard manner by the 
existence of a velocity notent ia^ $ such that 

““ u ■ ^ . 


It is therefore possible to 
reformulate (l) as a problem 
with elliptical type linear 
boundaries (figure 3) 


^ - 0 (fl) ; 

3 $ 

« - « < r > 


$ continuous 


continuous 





Note: As the boundary conditions are the Neumann type, it is pra- 
ctical to determine the potential at a defined point of the trail- 
ing edge 

$ » 0 (BF) 


1.2. Lifting Case 2-D 

An obstacle brought to light and the boundary of which may not 
be differentiated (point of reflection) can lift . 

The^lift (C z ) is introduced artificially by the Joukowski con- 
dition Ujjp • 0 (refer to GERMAIN (l)). In this case, (l) should be 
added to an additional scalar equation at the trailing edge 

J(u) - 0 (BF) 

This constraint makes it possible for a circulation S, to be in- 
duced around the obstacle, depending on the velocity at infinity and 
the shape of the body. 

A feasible modeling of the Joukowski condition imposes the 
equality of the pressures on both sides of the singular geometrical 
point (weakening of the condition u *g which is impossible to cal- 
culate numerically). By applying the law of Bernouilli 

P * f(|ui 2 ) ■ I - i 2 the condition of Joukowski is written 

I “—| 2 


l*l» 


i -*-|2 

,u 'EF 


(3) 




If (3) is added to ( 2 ), then a cut should be made (c) origina- 
ting at the trailing edge (BF) up to infinity (Figure 3) 

.Important Note : On figure 

BF 2 the second stop-point is 

not located at the trail- 
i edge : the fluid by- 
pi -es the obstacle, where- 
as 1 figure 3» the condi- 
tion of Joukowski necessi- 
tates that the fluid does 
not by-pass the obstacle. 

♦ ♦ 

On the other hand, if by starting at a point P gC the ob- 
stacle is by-passed and we return to point p“e C~ which is geome- 
trically mixed, then we have the relationship (4) 



<|>(P + ) 


<*>(P~)+£ , VPe(C) 


(M 


where & designates the unkown circulation 


iff! n ! jm-P 


Taking (3) and (4) into acoount t the formulation of the lifting 
problem analogous to ( 2 ) is written 


(5.1) - 0 (fl) • 

(5.2) $ + - l (C) 

(5.3) . | VI 2 - | V| 2 (BF) 

(5.4) - 8 <n r-r p ur„ 

(5.5) $ - 0 (BF) 


4 > discontinuous 


u continuous 


It may be noted that the non linearity cf ( 5 ) is due to the Jouk- 
owski condition and that the solution of the problem ( 5 ) is the 
couple where $ is a function and j, is scalar. "" 

In the tridimensional lifting case, « discontinuous sheet (ND) 
should be introduced at the beginning of the trailing edge line, fol- 
lowing the bisecting plane up to boundary and generated by the var- 
iations of the circulation in enlargement .( Figure 4). 

N 

•ND - u (ND) . 
j-1 J 
N 

LBF - J (BF) j 


Ip < bf J< 


Figure 4 


The coiling of this sheet for reasons of calculation time 
is left out and the formulation of tbo 3-D problem analogous to ( 5 ) is 
expressed : 



A6 - 0 


$ discontinuous on (ND) 


(A) { 


4> + (Pj) - 

l^(Qj) 1 2 


<r(p^)+*( yj ) 
- N(qT)| 2 
O') 


|u| continuous 

V Pj € (ND) j 
(LBF) 


(6) 


It may be noted that the solution 
(ift.fc) vhere 0 is a function and £ is a 



of problem (6) is the couple 
function of the erlargment. 
Finally, the formulation ( 5 ) 
is generalized in the case of 
a lifting flow around a multi- 
body (MC), by introducing K 
cuts originating at the trail- 
ing edges of the K— bodies up to 
infinity as is shown on fig- 
ure 5 (Example of a hyper lift” 
ing force (leading edge + pri- 
mary + flap)). 


Figure 5 


The problem at the boundaries to be solved is then: 


Find £.0 \ . _ 

**2 * *3 » * • • / solution of 


Ad> ■On** !-► 1 

‘ » r discontinuous j u I continuous 

♦ (p+ )-$(P ~)H V P £ c. i-1 2 3 

l^(Q + )| 2 

& e 6 (H 


(5) 


M 




VQe(BF). i-I ,2,3 


2. - FEASIBLE MODELING OF A SUBSONIC COMPRESSIBLE IDEALIZED FLUID 


2.1. 2-D Non Lifting Case 

As the flow is assumed to be irrotational , the compressibility 
model i? the isentropic type (refer to LANDAU-LIPSCHITZ (2) and the 
flow is controlled by the equation and the boundary conditions ( 7 ) 



p - p 0 <i 

7au ■ 0 


; u continuous 


ni 

Y+1 



1/(Y-1) 

) 


(0) 


pu*n " g 


(r> 


( 7 ) 


where p designates the fluid density 

y'the ratio of specific heats (y»l ,4 in the atmosphere) 
the critical velocity 

so that if we set 


P 


o 


k 




1 

Y-1 


then the law of compressibility is written 


P * O-kjuj 2 ) a 


If we introduce the velocity 
it is possible to reformulate ( 7 ) 
EAR boundary problem (8) 


potential $ by using V a u ■ 0, 

as a quasi-elliptical type NON LIN- 


(8.1) 7*p V(f> • 0 ; $ - * 

| (8.2) p - (l-klV*| 2 )° 
icS.S) p^-g 
*1 r 2 " 0 


u - * ♦continuous 


(«) 

< r i> 

(r 2 ) ; r - r, u r 2 


( 6 ) 


(8.4) 


In the compressible case. It Is interesting to add the local Mach 
number given by 


M 2 


2 r i V4>i 

ytl Li- |v«| 


2 



Furthermore, in the subsonic case, we have at every point of the 
fluid occupying the domain, the relationship ^2 < ^ 

2,2. 2-D Lifting Case 


Extension to the compressible lifting case does not present anv 
particular problem with respect to the incompressible fluid. The for- 
mulation is given directly by 


(pV0) ■ 0 ; $ discontinuous 

p - a-k|fy|y f!J) 

,-r (C) 

- i?n 2 (Bn 


JL continuous 


( 10 ) 


3. - MODELING OF THE POTENTIAL TRANSONIC FLOW FROM A COMPRESSIBLE 
IDEALIZED FLUID 


IXL 


3.1. Equations 

A characteristic of transonic flows is in the presence of shocks . 
The condition of irrotation necessitates that their intensity is small* 
^2 <J 5 where M is given by (9). 

If the local Mach variation is observed in a transonic case, it 
may be seen that there are what is called supersonic zones where the 
local Mach is higher than 1 (m > I) and subsonic zones where the local 



Figure 6 


In a transonic state through shock, the flow must satisfy the 

RANKINE-HUGONIOT conditions ?2) • .* .* + 

' ' » _ ^ «► — 

CP u *n] * [pu*n] 



with 


u *s continuous 

+ the region after the shock, - region before shock 
s » unit vector of flow direction 


n, 


orthogonal at * or in the direct 


sense of figure 6 -b. 

A characteristic of the transonic flow is that the fluid velocity 
may be locally discontinuous when passing through a shock. 


Let us now consider the equations and boundary conditions of a 
transonic compressible fluid. They are the same as those for a sub- 
sonic compressible fluid. 


V»pu ■ 0 Vau ■ 0 ; 
p - (l-k|u| 2 ) a (fl) 

[pu*n] + ■ Cpu'n] 

p • « <n 

The introduction of <p , however, by using Vau ■ 0 leads to a mixed 
elliptical type non linear ( 13 ) boundary problem * -hyperbolic 


03.1) 

V»pV$ . 0 0 

continuous 

03.2) 

P - (l-k|V*| 2 )° 

(«) 

03.3) 

[p^<j)« n ] + - [pfy. 

’nJ 

03.4) 

3d> 

p ^ ■ g (0 

ij discontinuous 


Fundamental Remarks 


u may be discontinuous 


( 12 ) 


1) There is no unlquity theorem for the solution of ( 1 3 ) in the 
transonic case, 

2) The conditions of discontinuity through the shock shall be im- 
plicitely satisfied in the variational formulation of ( 13 , 1 ), 


Figures 7» 8 show two possible Mach solutions in the case of a 
flow around a circle for „ 45 




DECOMPRESSION COMPRESSION 

SHOCK SHOCK 


PHYSICAL 

SHOCK 



Figure 7 


Figure 8 


The equation (l3»l)» where P ia given by ( 1 3 • 2 ) , is elliptical 
(hyperbolic resp.) in the regions of fl where the flow is subsonic (M < 1) 
(supersonic resp. (K>I)). The solution of figure 7 contains 2 shocks 
whereas the one on figure 8 only contains one. The latter is physi- 
cally acceptable, whereas the solution with a double shock, including 
a decompress. shock, violates the laws of thermodynamics. (refer to LAND- 
DAU-LIPCHITZ (2)). 

Accordingly, the formulation (12) or ( 1 3 ) is physically inade- 
quate. In order to prevent the appearance of non physical shocks, a 
condition of entropy must be added to 13* In the methods of finite 
differences, the condition of entropy is satisfied by introducing 
(M>l)of decentered differences or an artificial viscosity (see MURMAN- 
COLE (3), JAMESON ( 5 ), BAUER-GARABEDIAN-KORN (4) into the supersonic 
zone (M> 1) . 

In the finite elements techniques, the condition of entropy is 
treated as an added constraint to ( 1 3 ) or by a technique of artificial 
viscosity, similar to the finite differences, by modifying the equation 
locally in the? supersonic zone (l3.l). 

3.2. Tho condition of entropy formulated as a constraint 


During the passage of a shock wave, entropy increases and we show 


that in the case of a potential flow, this condition may be translated 
by a decrease in velocity through the transonic shocks. This charact- 
eristic applied to a monodimensional flow is translated by 


- u < 0 


(14) 


where u designates the velocity after shock, 
u" designates the velocity before shock, 
(Figure 9) 


If u ■ then (l4) 



d 4> < . , 

dx* 


(15) 


Figure 9 


By analogy ( 1 5 ) in bi and tridimen- 
sional becomes 

< +09 or < K with constant 0^) 
K to be selected 


In the variational formulation of the transonic problem, we shall 
consider a small shape of (l 6) given in (17) obtained by integrating 
( 1 6) by parts. 


- f $<fr#u>dfl s K fwdfl ¥ (o e 
J a Ja 

where 3 (fl) ■ {w|we (fl) ; wiO} 

& (fl) ■ {u(b)eC (fl) , supp (J compact) 


(17) 


It is important to note that in ( 1 7 ) only the derivates of first 
order, more accessible in a finite elements approach, are shown. 


The transonic formulation selected in this case is: 


^•p $$ ■ 0 

.P - (l-k|fy>| 2 ) a 
[p$4>*n] + ■ [pfy>»n!l 
p £" 8 
A<|>< K 


(18) 






3.3 The Condition of Entropy Formulated by the Artificial Viacoslt 


In roferonce to M.O. BRISTEAU (6) and to JAMESON ( 5 ), ( 1 3, 1 ) may 
be rewritten in 03)^ in a local reference marked (n/s) or j n | is 

the unit vector of the flow direction and the — ~ » s 

perpendicular orientated in the standard direBtion on ® ■*" 

figure 10. , 


in 


Figure 10 


p * _e_- (1 -» 2 ) i\ . 0 

3n 1-kV 3s 


2 . 2 
u +v 


v 3 u 9 

V 3x V 3y 


+ Z . (m . (2- £_ ) 
V 3y ’ ' v; A '3s * 3n ; 


In this form, the elliptical or hyperbolic characteristic of the 
equation appears, depending on whether or not V is smaller or larger 
than 1. In a similar manner to the decentering practiced in the fin- 
ite differences, the operator of artificial viscosity is added to 
(13.1) 

E(W - - ^ 0 > <„ 

with (|u|^-l) + ■ sup (0,|u| 2 -!) and the transonic formulation selected 
in this case is given by 


-V«(p7$) + vE(4>) - 0 ; v > 0 
P - (l-k|V$j 2 ) a 





Note t v parameter of viscosity >0 depends on step h of the triangu- 
lation of the domain in numerical applications* 

Other artificial viscosity operators mentioned in (5)» (37) es E 
in (21) have been tested numerically and give very close solutions 


E(*) - - «|u| 2 -1) + |A$|A*) ... 


( 21 ) 


3.4 Lifting case 


Extension to the transonic lifting case does not present any 
problem with respect to the compressible subsonic fluid* The formula- 
tion is given in (22) from (18) 


V»pV$ ■ 0 ; $ discontin * ; u discortin * 

p - (l-k|V$i 2 ) a 


CpV$»n3 + ■ [pV$»n] 

(ft) 

A<*>< K 

4> + ■ + A 

(C) 

1V1 2 » 1V1 2 

(BF) 

TS m 8 

(D 

■ 

0 

(BF) 


( 22 ) 


4 . - MODELING OF AN UNSTEADY INCOMPRESSIBLE VISCUOUS FLUID 


If ft and T designate respectively the domain occupied by the 
fluid and its boundary* the latter obeys the Navier' Stokes equations 
without dimensions* increased by boundary and initial conditions* i*e* 


/2? 


14 



(A) 


(23 


( 23 . 1 ) -y- - vAu + (£•?)£+ V p - 0 

( 23 . 2 ) 5 «u . 0 

( 23 . 3 ) u - 2 

( 23 . 4 ) u(x, 0 ) - u° 


where u is tie fluid velocity 
p is the pressure 

V Is the fluid viscosity ( v»i /Re with Re » Reynolds number) 
and « 0 specified; ^«o if p.p represent a wall (condition of 
2 adherence P 

If r«r represent "infinity" 

* njj oo 

An example of external flow around an airfoil is given on figure 

11 . 



Figure 1 1 

In the steady case (23) is reduced to 


— vAu + ( U *f)u ♦ ■ 0 



0‘; 

V*u - 0 


u ■ 2 

(0 


( 23 ), 


In the unsteady case, tho fluid is controlled by a system of 
equations wit 1 nnrnboli c typo no n linear partial derivatives, where- 
as in the steady case, it obeys a system of equations with elliptical 
typo non linear partial derivatives. In 


I 


Z] 
i 5 


ZE2 


(23) and (23) s the main numerical difficulties ire the condition of 
incompresslblity* the Reynolds and the non linear convection* 

5. - THE STANDARD ITERATIVE METHODS FOR SOLVING EQUATIONS WITH QUASI- 
ELLIPTICAL NON LINEAR PARTIAL DERIVATIVES 

5*1* The Model Problem 

For reasons of simplicity* ve are interested in the solution of 
the non linear Dirichlet problem (24) 


- A$ - T(0) - 0 

(fl) 

$ - 0 

(T) - Oft) 


with T non linear operator and (ft) a boundary of Rf 
5.2, The Fixed Point Methods 

The simplest algorithm to solve (24) is 


n*0 ; $ 


fiven 
to that 


r 


, n+I 


For ni 0 compute^” knowing by solving 




n+l 


T($") 


. n+I 

$ - 0 


(Si) 

<r> 


(24) 


(25) 

(26) 
(27) 


( 25 ) ( 26 ) ( 27 ) i« a converging algorithm for subsonic compressible 
flows (refer to GELDER ( 7 )* NORRY-DEVRIES (8), PERIAUX (9)). 

The Gelder algorithm in the case of a 2-D pipe. In this case, (qj 
is represented by figure 12. 



G 


n 


r» 


2-D Diver- 
ging Pipe 



In this case (r) ■ (PuD and (8) Is expressed 

• P 


(28.1) 

7»p 7 $ • 0 

(R) 

(28.2) 

p - (l-k|$ 2 $|) a 

(R) 

(28.3) 

f£-° 

( v 

(28.4) 

4>L - h(x) 

1 81 

<>V 


( 28 ) 


by 


The variational formulation ( 29 ) is obtained by multiplying (28.1 ) 
a test function 1 „ and by integrating by sections where 

U(H (R) 


H* (R) - {w«L 2 (R)|7wcL 2 (R) 


(29) 


f (l-k$ 2 $)°%»$n dx - 0, VutB 1 (0) , »L ■ 0 , • h 

'r *■ 

Let us introduce the functional ( 30 ) G q ($) and the space 

H 1 (R) - {*« H 1 CS1> | r ■ 0} 

os *s 

V« --W&T) d « 

Let us calculate f in the meaning of G&teaux (refer to VAINBERG 
(lO) ) f the derivative of Go at a point of hJ 8 


(30) 


Jij 6 0 <^ Xa ^*) • <*!£(♦*>,«♦*> 


The steady state of G 0 in is expressed in (31 ) 




i By approximating (33) from (29) t all tha steady points $ in 
o«(^ $*-hJp ■ 0 are •olutions of (28)* 

Wo can now prove the*uniquity of ( 34 ) 


min G (A) 


(3 1 *) 


in the case of the subsonic state by demonstrating that in this parti- 
cular ^se G 0 is convex * We have only to calculate for that in (33) 


C" - lit, 1- <G ($+X6$)) 

0 X-0 dX 2 0 

C"<$) - -f (l-k|$$| 2 ) a {fo<$«4> (^.?60) 2 } dx 

° (l-k|V$n (35) 


By refering to ( 9 ) M 2 - 2ko ( 1 -k | V4» | 2 ) “ l | V<*>| 2 wit" M local mach 


w* \ ✓ / ** » 1 T I . 

and using the identity j+.£| . |*||b| cos 9, ( 35/ is ytTxx ' r ' ' 36 ' 




i 


P(l-M 2 cos 2 0) |V6$| 2 dx 


(36) 


It is now easy to verify that if m< I in q -subsonic case) GJ 
is convex and that 


♦ - Arg min G (♦) 




os 


is the solution of (28) 


whereas if M>1 in 0 (transonic case) GJ is no longer convex and that 
is only a saddle point of G # , 


(37) 


A fixed point algorithm or quasi linearization is described in 
i) n*0 $° initialise en -A$° ■ 0 $°-h|p ■ 0 


the continuation {$ } is constructed by solving for $ n+1 e H 1 (J!) 

n n£ 1 



(37/ 



i 


The convergence of (37) in the hardest subsonic cases is obtained in 
10 maximum iterations . 

Note : The fixed point methods (37) are related to the gradient 
methods. In fact, (27) is a special case of (38) with p» 1 


$° given 

-A $ n+,/2 - T(4» n ) in £2 

/2 = 0 on r 

<j> n+ i * $ n + p(4. n+,/2 -4> n ) 


(38) 


but (38) by eliminating $> n+l /2 is expressed (38)' 


<j> n+l - $ n - p(-A)" 1 (-A<J> n -T(4» n )) 


(38) 


(38)' is a gradient method if T is the derivative a functional. An 
example of [j&)' described in (4l), consists of minimizing the funct- 
ional G 0 by a gradient method in the metric adapted to standard Hqs 

written out 11 si where designates the scalar product 

of the space 11 11 1 ’ «f f » = I 

of Sobolev Hl s . L 7f S* Vf 2 

Generally, if ^ 6 L 2 (fl) . Q , ^ defined by (39) is a dual ele- 


ment of 
L 2 (H)' - l 2 (Q) 


<G'($),5<|» 




pV<{i76<j> df2 . 


Let us introduce g e H 


os 


Then (4l ) consists of 
4 l .1 . 0 initialized in 


solution of (40) 
j $g.$6<j> dft - pV<f>- ^6<f)d0 , V6<J> c H 


1 

os 


( 39 ) 


(40) 


A(f» 


0 . <J>°-h| r - 0 


iZL 


we set h° = g° 


8° €H os (fi) calculated -Ag° - G * (4>°) 


(41) 


41 .2 . n -l knowing and g n e h* {<{> n+ *} Cg n+ *}eH* is construct- 

ed in two phases : 


os 


os 


Phase_l ! Calculate A * . arg G W n -X h ”) 

set - A" ° 

Phase 2 : Construct the new direction of descent 
To accomplish that solve 


Solve 


Calculate 


ans set 



(^♦l) is the version of the POLAK—RIBIERE conjugate (l) in the metric 

adapted H os . The rapidity of the convergence of (4l) is comparable 
to the one of the fixed point described in ( 37 ). 

The gradient method with metric adapted to the preconditioning 
shall play a fundamental role in the case of transonic fluids. 

5» 3» The Newton Methods 

Assuming this time T differentiable (24) may be solved by the 
algorithm (42) 


(42.1) 

<f>° given 


(42.2) 

nil , 4> n+l is calculated from. <J> n by 

-A4> n+1 -T'(4> n )4> n+1 - T(4. n )-T'(<J> n )*<J> n in « 

(42) 

(42.3) 

<£ n+ * * 0 on ' F 


a special case (p=i) of (43) used for the hardest cases 


(43.1) 

<t>° given 


(43.2) 

n2 j t ({) n+1 is calculated from by 



1/2 - T*($ n ).* n+I ^ 2 - T(4> n )-T'(4> n )*4> n in n 

^n+1/2 , 0 on T 

( 43 ) 

(43.3) 

4, n+1 . <{> n +p(<j) n+1/2 -4> n ) , P>0 



The treatment of the Joukowski condition, differentiable non 
linear constraint, is an application of (42). 

Example of the Floy Around an Airfoil 


(5) may be solved directly by using the discontinuous velocity 
potential in the form of an iterative procedure (Al) including the 
Joukowski. condition, the technique of decomposition 


of the potential (A2) (refer to NORRIE-DEVRIES (8)). 



A1 - If (c) designates an arbitrary lifting cut of the trailing 
edge 7 b f ) up to infinity (I^ , the velocity potential ^ Is a discon- 
tinuous function along (c). 

Let us introduce q ■ H* (H) designates the standard 

Sobolev space c 

H 1 («) - (v € L 2 (n) ; $veL 2 (fl)} 


then it is possible to give a variational formulation of ( 5 ) in 
/jj y tinder the space of H * (q ) defined in (44), 

£ c c 

H,[(n c ) - ( ve H 1 (fl c > | v| (BF) - 0 ; v| + -v| _ - . ( 44 ) 

c c 

Account taken of the continuity of the velocities along: (C) which 
is expressed in (45) 


-*•+ -►+ 
u *n 



9<J> I dd> I 

*=V'*H C ' 


(<.5) 


By multiplying 5*1 by a test function and by integrating in parts 
by taking into account (45)-(5«2)» the equation (5»0 is written in 
the variational form (4.6) 


L'V$*Vw dx - 0 Vu)€ H. 1 (fl ) , <pe H l (Q ) 

Jq 1 c c 

c 


(46) 


Assuming JK(£) 
is expressed in (47) 



the Joukowski condition 



(47) 


The algorithm (42) applied to this example is expressed in (48) 


1 i) Assuming £_° given ; $ is solution of (46) in H q (J 2) 
ii) For ns 0 f<£ n ,i. n } being known, {Jt n+I ^ n+ l} ^ 

are calculated by the aquation 

* n+I - £ n - JK*' , a n ) # JK(£ n ) (^ 8 ) 

with JK'» 2(i>. + 7 6$ + - V?6<T), ^n+l solution of (46) 
in l 

iii) stop test on ^ satisfied, otherwise nsn+1, go to ii). 
The convergence is ensured in several iterations (5 maximum), 

A2 - The velocity potential $ is the linear combination ( 50 ) of 
two potentials 4^ and continuous potential and discontin- 

uous potential, solution of (49)np and (49) D 


NP = g 


<*»)> 


A0 r - 0 

(«) 

d<P 

lK Rm 0 

(n 

*rL + ’ *rL 

* i 


- &*i 


I (C) 

- o 


C<9), 


* ' ♦» * *♦« 


If l is selected so that (51) occurs 




( 51 ) 


JK(JL) - m\ 2 ■ - (VJ|2 , 0 

BF + BF" 

it is then easy to verify that (<f>»£) solution of (49) » ( 50 )» (51 ) is 
the solution of (5)* 

Assuming “ (w e H ■ 0} and H* (fy is the sub- 

set of h 1^q j of the verifying functions on the cut (C) 
c 

“I * - “I _ - 1. 

c c 


Then the variational formulation of the equation '49),. m is expressed 
in (52) NP 


whereas 


f Q dx - j gw dr vu) C n l Bp 
* np € h - 

the one of equation ( ^9 )r is given in (53) 

j~ %-Vuto.o v ^ h ' to , 

c c 

$R enjoy 


( 5 : 


( 5 : 


The solution of ( 5 ) is then given by algorithm (54) 

i) <(> m and ^R solutions of ( 52 ) (53) ? ^o initializes ;$° = 4> Np U°< 

ii) For 0 J (i ,$ n } being known {£ n+ *,0 n+ *} is calculated 
by the equation 


£ n+l = £ n - JK'' 1 (Z n )JK(i n ) 


(54 


A 




with 


JKtt n ) - |$$ n | 2 - | fo n | 2 

BF BF" 

jK*a n ) - 2{fa n .fy I + -^ n .^j > 

BF * BF~ 

mX* I . , A r>f 1 
♦ - ♦ m + » »„ 


iii) If a stop test is not satisfied by n+1 J nsn+1 
we return to ii * 

The convergence of (54) is ensured in 3 or 4 iterations. 

Remarks : (48) and (54) are generalized in the two following 
directions : 


Uoo 


-compressible subsonic and transonic equations 
-complex geometries : 2-D multi-bodies and 3-D airfoils. 



Example : Expansion of (54) 
to a multi-body in subsonic 
state; The domain ($) is 
shown on figure 13* 

If P is eliminated in 
8,1 by using 8,2, the pro- 
blem with boundaries to be 
solved is given in (55) 


Figure 13 


(55.1) 

( *x' a2)l> xx ^ ♦ ^xVxy ‘ 0 


(55.2) 

4>l + " $\ _ ■ 

(55) 


cT c. 

1 1 



-► 2 , , 2 


(55.3) 

jku) = j V 4 -! _ = o i=i ,3 



BF . BF . 

l x 

\ (55,4) = g on r , = 1^ u ( u p.) , <f = 0 *n !% = { BF 3 } 


5 - 


(55«0 may be reformulated In (55»*0 in the form (55,5) 


-A<J> + T(<J>) - 0 


(55.5J 


,itl»T(*> 


fy + $ 2 $ + 2<J> <P 4 

xxx yyy xy? 


a 2 ■ A+B I V 4 > | ‘ 


A - a 2 + j (Y-l)'|uJ 2 B - j (1-Y) 


The method of quasi-linearization developed in ( 37 ) is used in 
(56) to construct a continuation ($ n }€H*(ft) verifying 


Ww dx - g (i) dx + I T($ n ) u dx ■ 0 


* n+1 c^n) 

Voj«H ^ 2 - {ui c H* (H) | oj| p ■ Q) 

If $ n+ * designates the discontinuous potential of the velocities 
and (Z.) n the circulations around bodies (P ) with iteration n f n +1 
is expressed in ( 57 ) „+i I . ^ 

* n+l - C ♦ x.h *U (57; 

" 1“1 

where < t, n+I and <*, are solutions of (58) ( 59 ) 
y NP y Ri 


Vd> n+ ] 

y NP 


7wdx -f g <0 dx + j T($ n )uj dx. - 0 
J r Jn 


C leH C 2 ( n > ; WoeH^fll) j 

^r! V “ dx « 0 Vw € hJ (Q c ) 

♦R i eH i.^ c . ) iH » 3 

The expansion of ( 5 ^) is the shown by algorithm (60). 

(<J*p ) solutions of (59) 

1 i-1,3 

(<*>^p) solution of (58) with T E 0 

solution of jk (£°) = 0 
° 3 

< * > ° initialized $ 0 = ^ J 

i 1 



5*4* The Pseudo-unsteady Methods 

They consist of associating to problem (24) a problem depending 
on time (6l ) 


5 * - - T<(J>) - 0 

i 

$l r “ 0 

^(x f 0 ) - 4 > (x) 

0 


<P> 

(T) 


( 61 ) 


The solution of (6l) : boundary <p(c,x) is obtained by using a spatial 
approximation, substituting for a system of normal differential 

equations integrated numerically on interval (0 ,t( , T large* 

In the case of an undlfficult problem, an explicite scheme de- 
scribed in ( 62 ) is adequate to integrate numerically (6l) 


i) 

ii) 


<D° - * 0 

n 2 0 


.n+l_x n _ n 

$ d- - A 6 n - T(4» ) 

At 

.n+1 n 
4> "0. 


0 


( 0 ) 

<r i 


( 62 ) 


Examples of the Pseudo-unsteady Approach t 

1, - Solution of the Navier-Stokes equations (refer to FORTIN (l4)) by 
the Arrow- Hurwi z algorithm 


The variational formulation of ( 23 ) is given in ( 63 ). 


Va(u,v) + b(u,u,v) - (f ,v) V v£ J , 

A “ 


(63) 


f £ L («) 


where J Q - (ve C v - 0} (for simplicity we have assumed) 

witha(u.v) - 7u*7v dx r 


0) 


b(u,u,v) *A(u*7)u*v dx 
(f » v) « /^f* v dx 


N, dimension 


of 


space 


Let us .now consider the discrete problem ( 63 ) . associated with 


V < P /’ V (P 2 )N « P h £P l and 


c| h € P 1 wliero designates 6ho Poly- 
gons with degree k 


^vv^Vs vH> . a ^ 

< H-V • 0 


(63). 


The Arrow-Hurwicz discrete algorithm substituted for (62) d may be 
described in (64) in the form of an explicite scheme. 


i) 

U) 



u^,p^ initialized 

nil known > (u^Jis calculated in (64,1) 



V w h e P 2 ; K ■ At 



K(f,w h )(64.1) 


/ * * * \ . f n* * 1 ^ j 

(in) n* I , lp h } & (u^ } known # P h is computed in <64*2) 

■ * v pl 

(iv) Convergence test on (u? +I ,p n+I ) not satisfied, do 
n=n+l, go to ii), « h 


(64 . 2 ) 


Note s The explicit numerical scheme described in ( 63 ) is rela- 
tively easy to program and economical to place in the computer core. 
Nevertheless, the conditions of stability connecting ^ . K and h has an 
industrial constraint. Furthermore, the numerical simulation of separ- 
ated flows, relatively hard case, requires several hundreds of itera- 
tions. 


2 ..Solu tion of Potential of Small Perturbations in Transonic State by 
Hnite Differences (Refer to I. A. ESSER ( 13 )). 

To the non linear system (63*) 


Fj » 7au ; u ■ (u,v) 

F - a u ♦ v ; a -0-M 2 - (Y+I)M w u) 
2 * x *> 


( 63 )' 






we relate the hyperbolic system (63") with suitably chosen boundary 
coalitions. 



F 


1 


3 v 

57 


bF. 


b > 0 


(63)* 


where b equal to the initial unity in H, YOSHIHARA (12) is optimized 
by taking b - j ct| to accelerate the convergence velocity of an ex- 
plicit scheme of second order of the Lax-Wendroff type. 

In the case of a very "hard" problem, it is better to use an 
implicit integration scheme to solve (6l) described in algorithm (65) 


i) “ *o 

ii) niO 

4 > n *|^> n _ - T($ n+1 ) - 0 (<>) 

St 

$ n+l . 0 ^ 


(65) 


at each step At a non linear (66) type ( 24 ), but better conditioned, 
non linear problem must be solved. 


/Id *\,4. n+ I 

<37 " 

- T($ n+1 ) - 

At 

(SJ) 

(66) 

A n+! - n 

♦ r 0 


<n 



6. - THE FUNCTIONAL LEAST SQUARES METHODS 

6 . 1 . Relationships between a Least Squares Method and an Optimal Con - 
trol Problem 

A least squares type formulation related to a model problem ( 24 ) 
is given in (67) 


min I |Av + T(v)]^ dft ■ min || Av+T(v) || ^ 
V€ V 'S2 V«v 2 

Ilf II 2 - f l.f| 2 

L h 

o , 

and V a functional space V" (ft) for example,, 


( 67 ) 


Assuming £ now the solution of the boundary problem ( 68 ) 


-AC - T(v) 

Sip ■ o 


( 68 ) 


( 67 ) is then equivalent to ( 69 ) 


min 

v«V 


f |A(v-C 


■o I dfl 


with £ « 5(0 via ( 68 ). 

By referring to J.L. LIONS (l5)» it is obvious that (68)-(69) 
has the structure of an optimal control problem where 


°0 v designates the CONTROL vector 
3) 5 designates the STATE vector 
Y) (63) is the STATE EQUATION 

gx the functional ( 69 ) is the function of cost or 
'criterion 

From ( 68 ) ( 69 ), it may be seen that other formulas are possible 
by selecting a different cost function. We may. for example, consider 
the optimal control problem (68) ( 70 ) 

Min I |v-£|* d ft 

vcV ■'ft (70) 

with 

5 I £(v) via (68). 

(68 ( 69 ) and (68) ( 70 ) shall give a solution identical to the 
solution of the model problem (24), but with different converging 
velocities . Furthermore, the choice of a least squares method is very 
important on the numerical level. In fact, a standard which is in- 
appropriate for the state equation ( 69 ) appearing in the cost function 
may lead to a slow convergence. A sound choice of the cost function 
with respect to non linear Dirichlet problems of the second order is 
discussed in paragraph 6.2. 

6.2, The Least Squares Method in a Particular Functional Space 

L: t us introduce in (70 (72) the Sobolev spaces required for the 
stu<2y of the model problem (24) 
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H*(ft) *{$e L 2 (ft) , V$eL 2 (!))} (71) 

H J( fi ) “ (Oc H'(ft), $| f . 0} 


( 72 ) 


(73 


'W 1 

‘ c H 1 (£ 1 ) 


+ [ 


$,$2^ + J V$j *74>2 dfi 


n 


Ik II 


H (Cl) 


f t 2 da ♦ [ (V4>! 2 an 

J n Jr 


Q 


(7'» 


H 1 ^) 


. s a sub-space of H (Q) . Consequently. If A is limi' 

ted, H* (Q) is a Hilbert space with scalar product (75; and correspon 1 
ding standard ( 76 ) 


( 4 >j *$2) t 
1 1 iT <rt) 


| 7$|*7$ 2 


Ik II 


,1 


( 


Assuming H ^ (ft) ■ (H^(H))’ the dual topolocial space of * 

By observing that 2 2 the inclusion (77) is permissible 

(l/ai)) 1 - l/(D) 

c L 2 (ft) c H -1 («) (77i 

Furthermore, the application ^ * $ 2 is an isomorphism of 8 n (H) 
in 1- <*,•> designates the bilinear shape of the duality 

between.,-1 and 1 defined in ( 78 ) by 

tr<n> 


h ;w) 


f 1^1 


2 dO) l/2 , 


(75 

(76 


.1 


h * <n> 


<f 


-1 


■ I 


f $ dx V f € L 2 (fl) ; V $€ H\f 2 ) 

a ° 


( 78 ] 


then the topology of H (ft) 
using ( 76 ) (78) 


is defined by ||.|| in (79) by 


<f,»> 


II £ II ■ fup 

4cJ<n)-(o) ||o|| H i (!)) 


(79) 


Refer to LXONS-MAGENES (l6) f NECAS ( 17 ) for more result s and 
characteristics relating to Sobolev spaces. 


By using (79) the best formulation in the direction of least 
squares for solving the model problem (2h) is given in (80) 


Min (J Av+T(v) || . j 


( 80 ) 



By introducing £ c H* (ft) solution of (68), then (80) takes the 
form (81) 0 


Min || A(v-£) || . 
vcH 1 ffl) 

o 


(81 ) 


By expanding ||M V -£)|| . into (82) 


<Mv-£),0> 


° o 

and by applying the Green formula to (82), we have 

|<A(v-£),0>|«|f |-<v-£>0 dr - f Vfv-£>$|> dft| - |f V(v-£)*$$ dfi| 

Jp3n jjj 

and (82) takes then the final form (84) |j $( v -£) «$<|>dfl| 

ll w >■* II _ _ ' n _ I! „_r'l 


( 82 ) 


IlMv-OlL, 


0 c Hi(P)-{0) 


lidl H l 


liv-Cil 


H o (f2) 


(83) 


(84) 


The least squares method in H** 1 ( 80 ) is, then, equal to an optim- 
al control problem , if,* 2 1 

“S$) (v> " W/ <v - £) i •iLi.V™ ■* * 1<v » :«5] 


|vcH;(fl) J C1 ' ' solution ' e « 

6,3, Iterative Solution of an Optimal Control Problem by a Conjugate 


radienl 


The Polak-Ribiere (ll) version of the conjugate gradient is used 
to solve ( 85 ), the algorithm of which is composed of 3 steps* 

i) Initialization 


0 I 

v 


given (for example solution of -Av° - 0, v°l ■ 0) 
gradient of J(v) in H 1 r 


is calculated in ( 87 ) 


We set t 


-Ag° - J’(v°) 

8°l r » 0 


(87) 


z # « g # 


( 88 ) 


, jThen for n - 0, assuming v n |6 n fJ! n known, calculate v n+1 , 

g n+ , z n+ by 

ii) descent 


Calculate X ‘ 
and set , , 


n+1 n ,* n 
v *= v ->.z 


(89) 

(90) 


iii )Cor»struction of the new descent direction 




Define g « ^(D) solution of problem (91) 


-Ag 

n+1 


n+1 


J'(v n+1 ) 


( 91 ) 


g 


then calculate 


4 > n+1 in (92) 


n+1 


and define z 11 *'*' in (93) 


J^y^ 4(g n+1 - g n ) d n 

~~T |vg n l 2 an 

In 


(92) 


'a 


n+1 n+1 . n+i n 

2 - g + y z 


(93) 


do n=n+l and go in ii) 


The two important points of the algorithm (86)-(93) are s 
l) - The problem of minimization to one variable (89) solved by dicho- 
tomy, or the Fibonnacci method (refer to "GOLDEN SEARCH" in P0- 
lak (ll)). 

Z) - The calculation of g n+1 from v n+1 requires the solution of two 
Dirichlet problems at each iteration (68) with v = v n+1 , and 
(91). 


The point ( 91 ) is detailed below. J* (v) is calculated in a 
standard way (derived from a functional in the meaning of Gateaux 
(refer to VAINBERG (lO)) in (9*0. 

Assuming 6v e H I <J’(v),fiv> = lim J(_ v . + t6v )- J (v) 

0 t+o c 


(94) is expressed by using ( 85 ) 


t^O 


<J' (v) ,6v> = 


V(v-£)?6(v-£)dfl 


8 


(94) 

(95) 


By differentiating ( 68 ) 6 £ e H 1 (Jl) satisfies ( 96 ) = t'(v)»5v 

6 C| p = 0 (96) 


Using (95) and ( 96 ) the final calculation of J'(v) is given in (97) 

< d* (v) ,6v> = J V(v-C)*V6v d8 - <T'(v)*6v,v-C> . (97) 

We recognize in (97) j> ( v \ £ H _1 linear functional defined on 

by (98) 


4> 


8 


7(v-£)*V<J> d8 - <T' (v) # <f>,v-£> 


(98) 


Then g 11 is the solution of the variational problem ( 99 ) 
n ,,1 


g e H‘(8) 


( 99 ) 


\ [ Vg n *^dfi = f V(v n -^ n )V4) d8 - 

)n 


<T'(v )<j),v -£ > , v<p £ H o (fi) 






with £ n solution of 


(soqx e n c H^<n> 


( 100 ) 


m 


The algorithm (86)-(93) shall be used systematically in the ap- 
plications of T to transonic flows and the the Navier-Stokes equations* 

7. - THE LEAST SQUARES METHOD IN H" 1 APPLIED TO TRANSONIC FLOWS 

7*1* Subsonic Non Lifting Case 

By retaking the problem with limits (8), the non linear operator 
T is given in (lOl) 




( 101 ) 


By retaking (85) with T in the form (lOl) the least squares for- 
mulation in H“ of (8) is given in (102), (l03) 


Min 

4>eV„ 


I f 

1 Jn 


dx 


( 102 ) 


g 

With V « (v£ H 1 (ft) | vi “ 0 } 

Irj 




dx = 


in 

V g « {0ev| p (4>) * g| r ^> 

C = 5(4>) via (103) 

P(4>)V<f>* Vw dx - I gudr 

ft J r„ 


Vli)£V. 


The physical interpretation of (103) is given in (104) 

AC «= V-p(4>)V(4>) i n ft SeV 


(103) 


(104) 


7*2. Transonic Non Lifting Case 


p(0) 


3n‘ T, 


u 


g " 0 * 


In the case of a transonic flow (l 8 ), in order to prevent non 

physical decompression shoc ks* a condition of entropy, which may be 
treated, must be added to ( 102 ) ( 103 ) 

-either as a linear constraint of inequality ( 105 ) 

A 4>< K . ( 105 ) 

In this case, a penalty functional of type (l06) must be added to 

+ .2 


(102) 


(Ap-K) 


dx where 


+ 

(A4>--K) - sup(0,A<fi-K) 


( 106 ) 




Ml 




leading to the least squares method (102) with penalty 

P 

V I 2dX + U I Q I ( ^~ K) + l 2<lx 
with y > 0 and solution of ( 103 ) 


(1° 2 ) P 


Two possible alternatives of (l02)p are given in (102) and 
(102) RO with K=0 i this is a least squares method with regularization 

Mini{[|V£| 2 dx + yf |(A<j>) + ! 2 dx ( 102 ) 

J Q J fi 1 } R1 

with M>0 and £ solution of ( 103 ) 

Mini} JVC! 2 dx. U| f M * 2 <!*■>„,[ [S.Sf 2 * 

4>eV 'ft ‘Jq 

with ( )+ = positive intensity of a discontinuity = sup (o.( )) 

-or by artificial viscosity . in this case the functional (102) re- 
mains unchanged, but £ = £(«{,) v i a ( 1 03) y 


V£* Vu dx = 

fi 


p(4>)V<}>* Vo) dx + v a (4>)uj dx - 
ft J Si 


Vw e V ; £ e V 



g Li dF 


U°3) v 


a defined in (19). 


In the applications, (l02)p is preferred to (102) due to the 
sensitivity of P in ( 102 ) . In both methods, U is adde§ to obtain a 
same magnitude for both terms of cost function. Finally, it is also 
possible to combine the regularization given in (102)^ with the arti- 
ficial viscosity (l 03 )y to eliminate the numerical instabilities in 
the region of shock. 


7. 3» Transonic Lifting Case 


By using the notations of 1.2 and by referring to the lifting 

flow shown on figure 3 » the circulation l of u * around an airfoil 

is in general ^ 0. Thus, t is discontinuous and a cut (c) (figure 3) 

must be made, * 

Assuming SI « ft - (C) & JK : R -*• R the function defined by 

JKOO = | ($$.) + | 2 - | <$§a) _| 2 (107) 

x BF BF 

where is the solution of the physical problem ( 107 ), (l08) 


% 

3n 'r 


in 


u • n 

00 


C + 


- * 


SI 

Mi , 
’ an'r 

p 

M. I 

9n+' 


= l 


+ M 1 
c + 9n c 



( 108 ) 


The method of decomoosition described in (49)-(5M is applied 
to the lifting transonic case. 

By selecting the least squares method with regularization (l02) R 
or artificial viscosity ( 103 ) » the following algorithm R 

is searched for in the form ( 109 ) 

“ ^NP + ( 109 ) 

where $jjp represents the NON LIFTING compressible part of the potential 
and the LIFTING. i.ncompressible part of the potential, '♦“r discon- 

tinuous on (c) is solution of (*»9)np# is solution of two iterative 
algorithms (TRANSONIC FLOW CHART 1-2 J * 

1, External fixed point algorithm defines £ solution of the non linear 
monodimensional equation (llO) 

JK(S.) - 0 (HO) 

2, Internal conjugate gradient algorithm gives 4>jjp» at Z fixed, as solu- 
tion of optimal control problem (ill) ( 112 ) 






INITIALIZATION (£,£”) (g ,h u ) 



TRANSONIC LIFT FLOW CHART 1 
ORGANIGRAMME TRANSSONIQUE PORTANT 1 








TRANSONIC LIFT FLOW CHART 2 
ORGANIGRAMME TRANS SONIQUE PORTANT 2 


3R 


limn nir-aiini—iHlil— Mr 










(Ill) 


.0+1 

Min r (<}> Np ) 

V v 


if i?Ei 2 d««J i«# t n *'> + r 

J n Jn 


dx 


VITH5 “ SOfryp) via ( ll2 )» * 4,n+I • ^NP +S,n+1 ^R 


♦where 


| dx * | p(^ + 1 )V^£ + »^ti> dx - 8 u dr V(i> e V (112) 

7.4. Conjugate Gradient Solution of* Non Lifting Transonic Problem 

For reasons of simplicity, we are limiting the problem to regulr 
arization (l02)p i.e. 


with 


Min J( 4 >) 


J( 4 >) - iff |V£| 2 dx + J | (A<|>) + | 

* J n Jn 


dx 


(113) 

(114) 


where V and £ defined in (102) (103). 


In this case, the conjugate gradient algorithm similar to the one 
given in (86)-(93) consists of three phases : 

1° ) Initialization 

<J>° is selected as solution of the incompressible flow i.e. ( 1 1 5 ) 

<J>° * 0 in ft 


< * > ° I r j " 0 • p lr'r 2 “ 8 


(116) 


of which the variational formulation is given in (117) 

JwdT Vojev ; <|P e v (117) 

r 2 ^ 

(if u is known on boundary ^ > P is also known by ^ * P 0 0“ k !u| ) ) 
since o° e V I s calculated as the solution of the variational equation 
(118) 


[ V'&» 


dx = <J ' (4>°) ,u» V 0 ) e V , g° £ V 


(118) 


Accordingly, we set h° 5 g °„ Now for 0,, assuming tp n ,g n h n as knov n . 

we compute ^n+I n+1 n+l i n two Phases* 

< >8 >h 


2°) Descent to calculate a^+I by minimizing the functional to one sin- 
gle variable (119) 


■ Arg min J($ n -Ah n ) 
A*0 


we can then set 


(119] 

(120] 


♦ n+1 - - X V 

3° Construction of the New Direction of Descent 

Define g n+, e V as the solution of the variational equation ( 1 21 ) 

$g n+ « 1 $u> dx - < J ' («f> n+ 1 ) ^i> Vide V , g n+, £ V 


'ft 


( 121 ) 


calculate the coefficient of conjugation n +j in the metric relating 

to v y ■ 


,n+I 


set then 

and return to ( 119 ) 


* n+l A, n+l n%, 

Q Vg »V(g -g )dx 

f l? 8 ”i “ 


( 122 ) 


dx 


n+l n+l n+l. n (123) 
h ■ g + Y h 


Note : Each iteration requires on the average 5 solutions of the 
Dirichlet problems : 


n+l 


-2 for the calculation of the gradient g in the good metric 
-3 on the average to calculate the optimal step 


Let us now expand the calculation of J'($ n+1 ) 

If <•,.> represents the duality between V* and V, by using (114) 

$£•$65 dx + vf (A<+>) + A64> dx (124) 

! - J ft /* 


<J' <♦>,«*> 


'ft 


a> " ® 

where §£ is the solution of the differentiated variational equation 


v lD£V, S 


(125) 


and 6p is expressed via the relationship 


P($) - ( 1 -k j | ) 




6p(<J>) - -2ha(l-k|$^J 2 ) Qt " l $4>.$5(|) 


(126) 


in (124) may then be expressed as a function of 6$ with (l 26 ) 
and ( 125 ) written with w ■ £ 

f dx - f p($)$?»fa$dx- 2 k f ( 127 ) 

J n 'n 

With (127) <J' (<{>), w> may then be identified with the linear 
functional 

u) ■* f pfo)$£.$w dx - 2ka f (p( 4 >)) , " ,/a (^<l>‘^ 0 (^*^w)dx + y[ 'A0) + Au) dx (128) 

Jn Jq J a 

From (121) (128) we obtain g n+1 from ^ n+I by solving 


f $ g n *.% dx - f[p(p n+, )^ n+ 1 *^- 2 ka (p(r l )) H/a (V +, ^ n+1 )x 

Jfl J Si ( 129 ) 

(fy n+1 •$*»)] dx + J (A$ n+1 ) + Aa) dx , ¥ wcV, g n+1 € V 

'n 

with £ n+I solution of ( 103 ) corresponding to <p - ^ n+I t 

8 . - THE LEAST SQUARES METHOD IN H " 1 APPLIED TO THE NAVIER-STOKES 
EQUATIONS 


8.1. The Steady Case 

8.1.1. Functional Least Squares Method of Steady Navier-Stokes Equa- 
tions 

In the following, we shall designate by ( 1 30) the scalar product 
<’. ,0 of two functions u,vc(H'(n)) N N atandln £ for the dimension of 
the space *' 


<u,v> - [ dx - l f $u.»$v. dx , V u,v € (H 1 (fl)) N 

; n i -iJfl 1 1 


(130) 


-*■ / \N ■+■ r iN 

“ " {u i’i-i - v ■ { Vi-r 

Let us define in (130 
2 


^ ■ {v € (H 1 (ft)) N , ^*v * 0 in D 1 v|j, ■ z } 


(131) 


Then the variational formulation of the unsteady ( 1 32 ) Navier- 
Stokes problem 


-vAu + (u*^)u + = 0 (fi) 

1-U - 0 (ft) 


(132) 


u| r * 2 


is given in ( 1 33 ) 


vf dx ♦ [ v» (u«$)u dx ■ 0 V vc Wj, ut Wj . 

in J n 

A least squares method of (132) ( 1 33) is given by the optimal 
control proble.i ( 1 3*0 


min (J(v) 




where "I in (134) is a function of y via the state equation ( 1 35 ) 

- a£+ ?tt — (v»^)v (h) 

H - 0 (R) ( 

1\ T - * 

of which the variational forsnulation is given in ( 136 ) 

vf dx ■ -f rj* (v*^) v dx: v n € w , wj 
J n ip 

It is essential to note that (135) ( 1 36 ) is a Stokes problem, 
acting as a pressure in ( 1 35 ) 

8,1.2. Conjugate Gradient Solution of (134) ( 136 ) 

The algorithm is composed of 3 phases: 
i) Initialization : 

Take for £0 the solution of the Stokes problem (137) 

-vAu° + $ p ° o 0 (ft) 

( 0 ) <1: 

*►0 ■+ 

U * z 


of which the variational formulation is given in (138) 


V dx * 0 V r\ e W , u° c K*- 

h 0 z 


(138) 


i 


Take for g «W the solution of the variational equation (139) 


[ d*-<J ' (u°),n> v « w , g° £ w 

Jq o 0 

and set ft 0 ■ g°. 

For n*0, assuming: as known, calculate 

ii) daacfiDi.aiiaafl (i4o) (144) 

X n ■ Arg min J ( u n -Xli n ) 

X>0 

7*' - u n - \«h n 


(139) 


by 


(140) 

(141) 


iii) phase of^cor a tract ing_the_nev_des cent _direct ion 

Take for g n+l e w the solution of the variational equation (142) 

(142) 


f r s **.% dx - <j'(u n *'),n> » » . ?♦'« w 

J n 0 0 

f $r'-W-?)dx 

n+1 _ 'fl 


Calculate n +t 

Y in (1 A3) 


Jjfg n | 2 dx 

The new direction of descent it”** is given in (l44) 


(143) 


h n+1 - l n+1 ♦ Y n+, h n 


(144) 

do nsn+1 and go in ii)* 

It may be observed that (139) (140) (l4l) are Stokes problems. 


8.1.3. Calculations of J* and of -*-n+l 

g 


By definition, the calculation of J* is given in ( 1 45 ) 


6j - <J'(v) ,<5v> - v 


^(v-£)*^(v-S) dx 


(145) 


'n 


It is possible to express ,as a function of <J V by using the 
differentiation of ( 136 ) 


<H6) 


v[ dx *•' n*R<$v*V)v + (v»V)6v] dx V n € W q ; 6? e W Q 

Since (v-0 t W Q f let us select n ■ v-£ * n 0^6) and 6v ■ n 
in (145) which is expressed t 


<J'(v),rj> 

vKcW o 



(v-$) *■ (n*^)v] dx 


047 ) 


To calculate g n+, t we must therefore begin by <j' (u r+l ) ,n> 

which requires the solution of the state equation ( 1 45 ) for -+n+I 
giving ^n+i (147) nay then be expressed in ( 148) y 


<j'(i n *').V - d, ♦ J[(t n * , -t n * l ).<i n * l .5 ) ; 

♦ (u n+l -5 n+I )* (n*^)u n+l ] dx 


(148) 


Finally g n+J is given by (142). 

In conclusion, each optimal control iteration requires several 
Stokes problems : 

• Stokes problem ( 136 ) to calculate the state £ n+ * from "u n+ * 

• Stokes problem (142) to calculate the gradient -*-n+! from "^ n+ l 

and £ n+ * ® 

•Stokes problem (l4o) to calculate \. 

Furthermore, an efficient Stokes algorithm shall prove to be a 
particularly important tool in the solution of the Navier-Stokes equa- 
tions vio the least squares method (134)-(136). Ite implementation 
shall be described later on in paragraph 8.3. 


5 1 


8,2. The Unsteady Case 

8.2.1. Formulation of the Unsteady Navier-Stokos Problem 

As was presented in paragraph 4, the unsteady Navler-Stokes pro- 
blem consists of (149) ( 150 ) ( 151 ) 


■|~ - vAu ♦ (u*^)u + $p • 0 

(A) 

(149) 

$-u ■ 0 

1 [ •+ + _ 

<fi) 

(150) 

u 1 p ■ z ; z»n dr ■ 0 

* r 

(O 


u(x,0) - u (x) 


(151) 


where the function z, given, may eventually depend on t. 

8.2.2, Quantification in Time of the Problem (149) (150) ( 151 ) 

Several schemes may be used to solve (149) (150) ( 151 ). For 
reasons of simplification, we are presenting two very simple ones 
with a constant quantification time step. 

8.2.2.1, Semi-implicit Scheme 

Assuming k « At the quantification time step. A semi-implicit 
scheme in time, which is very simple, is given by 


i) 


-+o 
u 


u c given 


(132) 


then for n 2 0 , u n+ *is obtained from by solving ( 1 53 ) 

1 -♦n A , 

u -u 

£ : " vAu 

( 0 ) 


^ n+l + tf 0 n+l 


-<u n -^)u n 


<n> 


- 0 


(153) 


-*n+ 1 1 l 

u | r - z • z((n+l)k) 


with u in ( 153 ) an approximation of u ' nlt ' where u is the solution of 
(149) (15l). It may be noted that in (153) ^ n+l is obtained from 'JJ*' 
by solving a linear problem , variant of the steady Navior-Stokes pro- 
blem 8.1 (here also the operutor S = -vA I* substituted by ^ Id 

Accordingly, it is necessary to develop an efficient k k ’ 

Stokes algorithm relating to in order to solve (l49) (l50) ( 151 ). 


Z2 



8.2. 2.2. Implicit Scheme 


of 


The simplest implicit scheme for solving (l^9) ( 1 5 3 ) consists 


i) 


■*o -► 

u “ u 0 given 


(15M 


Then for 


n i 0, 


•♦n+1 

u 


is obtained from by solving ( 1 5 5 ) 


ii) 


Ht n+ 1 '*' n , 

— _u - v^ +I ♦ 


k 

I 

V»|j 




♦ $p n+l - 0 (fl) 


yi+i 


o 

r -r +l 


(«) 


(155) 


-tn+t **T) 

It may ' «.* observed that in (155) u is obtained from u 
by solving a NON LINEAR problme, variant of the steady Navier-Stokos 
problem 8.1 (here also, the operator S « is substituted by 

g Id .It is from (155) that ve shall present a least 

5 k * "k squares method similar to that given in 8.1 for the 

steady problem. 


I 



8.2. 3. Abstract T-oast Squares Method from ( 1 55 ) 


In fact (155) 3 s a special case of a family of non linear pro- 
blems S Q (a> 0 ) (156) 

au - vAu + (£.$)u ♦ % . 1 (ft) 

* 0 


056) 


u|p ■ z avec 


( 0 ) 

f -*■ -¥■ 

= *n dr ■ 0 (D 
} T 


of which the variational formulation is given in ( 1 5 7 ) 


aj u»v dx + vj dx + i v.(u.^)u dx * [ dx, \*v e W ; u g w+ 

J n Jq Jq o z 


(157) 


By following 8.1 an optimal control least squares method of ( 1 56 ) 
( 157 ) is given in (I 58 ) 


z •• * n 


where t is a function of via the state equation (159) 
s v 


- vA^ + ^tt ■ "I -(v*$)v (ft) 

- 0 . (ft) 


t\ 


T 


(159) 


IT, acting as a pressure. 

8.2.4. Conjugate Gradient Solution of ( 1 58 ) ( 1 59 ) 

Tracing paragraph 8.1.2., the conjugate gradient algorithm for 
solving the least squares problem ( 1 58 ) ( 159 ) is given by 


i) u e W£ , given ( 160 ) 


calculate g° solution of the variational equation 


'msrjszm 


af g°*n dx + v[ $g°*$n dx ■ <J ’ (u°) ,n> , v n« w ; g° £ w 
Jn Jo o c 


( 161 ) 


and set ' h° - g°. 

Then, for 0, assuming u m > g ra ,h m i;s known, calculate Tm+1 ?hn+l by 

u > 8 » n 

ii) Descent Phase ,m . . , 

X - Arg mm J(u - Xh ) 

X>0 ^ 


•+-m+ 1 •♦•in , mf-m 

u - u - XTi 


iii) Phase of constructing the new direction of descent 

Define g n,+ ^ solution of the variational equation (164) 


( 162 ) 

(163) 


i ->-m+! ->■ , 

a g »n dx + 

if) 


dX " <J,( ’“ n+1), ^ > * V He W o , g n+1 £ w fl 


then 


af F)<lx+‘J ft*' •*(?*' -?)dx 

• . Jo__ in 

en af | g m | 2 dx + vf |^g m | 2 dx 

in in 


the new direction of descent ^-i+l is then 


h** 1 - | m+I + Y nH-l-gm 


do Rism+l and go in ii). 


(164) 


(165) 


The calculation of J'(u m+ *) is not detailed, as it is a trivial 
variant of 8, 1.3, 

In a similar manner as the algorithm ( 1 3T ) — 144), each iteration 
of (160) ( 165 ) requires the solution of several Stokes problems 
of type (149) without non linear term. 

.the Stokes problem S^ to obtain from u m+ * 

.the Stokes problem Sfc to obtain g®*! from + m+ i 

..the Stokes problems Sk *° calculate n u » 4 

A • 

\'ote : The algorithm -m (l6o)-(l65) permits the calculation 
n n+ l from -*-n as a result ' m( is initialized in ( 160 ) by -+n+l ,0 

u u V U ) U =■• ! 


8,3, A Rapid Stokes Algorithm (The Continuous Case), 

8,3.1, Summary 

Paragraphs 8,1 and 8,2 have demonstrated the necessity of devel- 
oping an efficient Stokes algorithm S q defined in ( 1 66 ) 

ctu - Au + $p - f (fi) (166) 

$.u » 0 (fl) 



for solving the steady and unsteady Navier Stokes equations. In 
(166) a=0 corresponds to the steady case, whereas a > 0 corresponds to 
the unsteady case. 

We shall show that the solution of (166) by following GLOWIN- 
SXI— PIRONNEAU ( 18) ( 19) (49) is reduced to the decomposition of the so- 
lution into a finite number of Dirichlet problems coupled with an in- 
tegral equation, 

8,3,2, Principles of the Method 

Let us note that by taking the divergence^ of the first equa- 
tion of (l66), we obtain an equation on pressure of type ( 1 67 ) 


Ap * 


(167) 


If we know ^ 

obtained by solving the 


then the solution ^ U »P) 
N+l Dirichlet' problems 


of ( 166 ) should be 
( 168 ) (169) 


Ap = f 
p! r - * 


(fl) 


au. - Au.= - 5-^- + f. (fl) i=l,. 
1 1 9 x. 1 ’ 

l 


( 168 ) 


,N N ■ dimension of the 

space (169) 


u i'r 


* z. 


Put wo don't know A ' 

The introduction of $ solution of ( 170 ) will make it possible to 
s i;t A , i„o. tho pressure trace on the edge, so that the constraint 
distributed u = o ^- 3 satisfied. 


~A<f> ■ u (jj) 

♦l r - 0 


( 170 ) 


In fact, by applying the laplacien A at ( 170 ), we obtain (171) via 

( 166 ) 


I 


-A(A<{)) - A($*u) » (Au) « Ap - 
2 

A 4> + aA$ « 0 (ft) 

♦l r “ 0 

34> 


(171) 


( 172 ) 


If we now select * so that 3 n * °* then after ( 172 ) and /-g 

fore > 0 . The application 34>Xi via ( 168 ) ( 169 ) ( 1 70 ) ^ “ 

lin-^^Pr 

|flr - ^ + b 


therefore v . u _ u . 

being affine, there is (A,b) (A lin- 3n I T ear operator, b constant) 
so that ( 173 ) occurs 


(173) 


Also, the (N+l Dirichlet problems ( 1 68 ) ( 1 69 ) coupled with the in- 
tegral equation ( 1 7*0 

AX + b - 6 (17*0 

give the solution (u,p) of the problem (l66). Let us point out that 
the good conditioning of the operator A is necessary to solve (17*0 
easily. 

8.3.3. 'Functional Support of the M*ethod 

To define (A,b) in ( 1 73 ) , it is necessary to introduce 


U ,/2 (D - (y e H 1 / 2 (D, lydr - 0 } 


• jv 

T 


(175) 


The method of decomposing the Stokes algorithm is then based on 
the following result : 

Theorem 8.3.3.1 . : Assuming XeH~ ,/,2 (D J assuming H _1 ^ 2 (r) h’^ 2 (D 

the linear operator defined by 


Ap^ » 0 (ft) ; 

p^ € H 1 (ft) p^ - 

»v H = _ ^ p x 

(fi) ; € (H 

= t *“x 

(«) ^ eH o (n) 

34> x 


A X " 3n 'r 



J 


(176) 

(177) 

( 178 ) 








Therefore A is an isomorphism of H '^(0 / R on H *^(T)and also the 
bilinear form a(*t«) defined by ( 179 ) 

a(A,p) - <AA,p> ( 1 79) 


/5 7 


1 / 2 

where designates the duality product between H (D and 

H - *^(r) is continuous, symmetrical and highly elliptical inH~^(r)/R 

The application of theorem 8.3«3»1» to the solution of the 
Stokes problem will now be possible thanks to theorem 8.3»3*2. 


Let (L 2 (J2» N ; & P 0 '“o’^o defined by 

Ap o -^f (fi) !P 0 el£(fl) 

au-Au « f-^p (ft) ; u-z e (H^ (ft)) N 

o o o o o 

-A(|) o - (ft) ; <J> Q c H^(ft) 


(180) 

081 ) 

082 ) 




‘Hieorem 8. 3. 3.2. : If (u,p) is the solution of the Stokes problem 
( 1 66) , then the trace A of p p is the solution of the lim 
variational equation (e) 


( Ae H' ,/ 2 (r)/K 


(E) 


(183) 


<AA ,y> 


<!!° ,y> VucH' l/ 2 (D/» 
3n 


The demonstration of these theorems is given in R. GLOWINSKI-O, 

PIRONNEAU (18). 

Notes : 

1) Theorems 8.3»3*1« - 8,3«3»2. show that the Stokes problem 

( 1 66 ) mav be decomposed into a finite number of Dirichlet problems | 

(-A) (resp. ald-A) (\+2 to obtain $ f N+1 to obtain r~*~ ^ when > is 

known plus the problem (E) ; 0 _ U ’^ 

3^ 

2) In the approximation 0 : h) or (E ) 3 n shall not occur expli- * 

citely duo to the Green formula applied in ( 1 S h ) if p is sufficient- 
ly steady. 


KSRffi,' 


* -~v 


» y> ■ J^v^c dx - j n d * 

r + ^ (184) 

■ ( $4> +u )*^v dx 
0 ° 

■» 

where li designates a steady rise of (j in (2. 

3) The main difficulty .lies in the fact that the operator A is 
not explicitely known. 


To overcome this difficulty! a new variational formulation of 
the Stokes problem shall be used in the approximation of (e ) t requ- 
iring a quantification into mixed finite elements. 

8.3*4. Mixed Variational Formulation of the Stokes Algorithm 

Let us introduce 


Wj - i{v,4>) € (H , (n)) N xH*(n), v| - z, f fy4)dx - [ ^.viodx, VWfiH 1 ^) 

h h 

W Q - {{v,4>}e (H^(n)) N+1 , | Vfc, dx m I V»V 03 dx Vu) € H* (ft) } fl 8s ^ 


It is easy to demonstrate proposition 8. *3. 4.1. j 


If {v,4>}e W*, then { solution of (186) 
z 


V (Q) ; <f> - !£■ - 0 (D 


(186) 


Ve have only to use the definition of W z and the Green formula 
in (185). Let us consider the variational problem (p) (I 87 ) 


Find {u,^} c W so that 
z 


a u.v dx+ ^u*^vdx=[ f*(v+$j>)dx, V{v,a} c W 
J Q J Q 0 


(187) (P) 


Theorem 8. 3.4.2. : (p) has only one solution ( u »4') where ip" 0 
and is the solution of the Stokes problem (166). The demonstration 
of this theorem is given in R. GLOWINSKI-O. PIRONNEAU (l9) shows that 
(?) is a mixed formulation which is interpreted bolow : 

. 2 1 
If V£ (H*(Q)? and T is sufficiently steady 3<J>£ H n H q (Q) an( j 

► -+■ I M so that the docompostion (l88) is the only ono. 

J ■ Va X € (h'(Q))‘ N 




( 188 ) 


In the formulation (p), insteady of setting directly V*v ■ 0 
wo try to set <t> = 8* which is equivalent in the continuous case, but 
not in the discrete case* The approximation of (p)h from (p) via 
the mixed finite elements shall be presented in paragraph 10, 


9. - APPROXIMATION BY THE FINITE ELEMENTS METHOD OF THE TRANSONIC 
FLOWS 


9, 1, Summary 

In this paragraph the approximations by the Lagrange finite 
elements of the transonic flows considered in paragraph 7 are briefly 
reviewed. Refer to the works of M.O, BRISTEAU (6), (20) ( 38 ) and R, 
GLOWINSKI and 0, PIRONNEAU (21 ) for more details. 


For reasons of simplicity, only external flows around airfoils 
shall be considered, 

9,2. 2-D Flows 

9*2,1, Case of Non Lifting Airfoils (profiles) 

The situation is summarized on figure 15* 




Figure 15 


The _transonic flow around a symmetrical airfoil P is without in- 
cidence (u^ parallel to the chord of the airfoil) and is modeled by 
the relationships of paragraphs 7*1» 7.2, The flow symmetry results 
in automatic satisfaction of the Joukowski condition (,L07). 



9*2.1, 1, Approximation of the Spaco V 


By retaking the definition ( 189 ) of the space V in 7*1 » 7*2 

V - {4> € H 1 (Q) nc°<n), $-0 at trailing edge ^ 

If d signates a polygonal approximation of the domain ^ occu- 
pied by the fluid and if is the set of triangles (Tk) or TRIANGU- 
LATION such that y in a standard vay 


“ T k V T j - ♦ 


(190) 


then V is approximated by the space of the finite dimension Vj, 

V h ^h eC \I T€P j c v Te Tf h , $.■() at trailing (l9l) 


'h» 

Similarly, if we define by (l9l)^ 

3<J> 


edge 


v h g - { »h tv J 


(191) 


e 


In (l9l) P designates the space of polynomials with two variables 
with degree sk. In practice, the numerical tests require k=l or 2, 

9.2.1.2, Approximation of the State Equation 

The state equation expressed in (103) is approached in ( 192 ) 


L VK “* ■ L P<V ** - f S„ % ir t 

h n r h 


( 192 ) 


V V hg ’ Vl V V h 


where g^ is a suitable approximation of g on the edge 

If k=l are Piece-wise constant over each T « con- 

quently, p/* \ is also constant and ( 192 ) may be calculated accura- 
tely. h 


If ks2 . are piece-wise linear and a numerical integra- 
tion of is necessary. We may proceed as follows 1 each 

is divided into 4 sub- triangles (Figure 16) 



On each Tj (jsl f 2 t 3 f i|) is substituted by a linear inter- 

polation P, . This approximation permits again an accurate integra- 
tion of (192) via the FORMAC system, A. LAPLACE (22), 

9. 2, 1.3. Approximation of the Cost Function and of the. Penalty Funct- 
ional 


We approach the cost function by (193) 

W ' T |„ 


(193) 


For the penalty functional, two approximations shall be consi- 
dered : k=l ; ke2. 

The linear constraint of inequality (105) is expressed in a weak 
form (19*0 

■Jo VK d * Or n dr 4 K L “h d * * “h ‘ »; (ISO) 

2h % 

where is the sub-unit of according to 

v h • { VV V 0) 


(195) 


If the bounded K is also defined with the weak meaning by the 
variational formulation (197) from ( 196 ) 


4 *oh- K <V 

♦oh ■ 0 "V 

H ° h . fr 1 
"35" *h|p <r 2h ) 

| ? “h dX ' 


-K| a), dx + 

J «h lr t 


h "h dr h 


( 196 ) 


(197) 


( 198 ) 


The constraint (19*0 is substituted by the discrete condition 


'In. dxs0 v “h eV h 


( 198 ) 


Lot bo a baso of Vjj produced by the functions of form N\ 


* ^ N 0is| w ithN T h . dim(v h ) 


(199) 


(M j ) - , V Mj e {nodtaoef Tf^} (Trailing edge noda) 


( 200 ) 



Figure 16 


We may then substitute for (198) the N. constraints of ine- 
quality (20l) 

q i • -L ’W'K *■ ° \ 


(201) 


One way of satisfying them is to add to criterion (193) the dig 
Crete penalty functional (202) 


P 1 23 " l l Q il 
icN, 


+ .2 


( 202 ) 


If k=2 , the base functions N. are not all positive, 
case may be decomposed as follows in ( 303 ) 


In this 




( 203 ) 


h *’h 


B 


1,2,3 


(N. , i c tops 1,2,3 of triangle T c 


tt> ** * 5, 6 _ 
h 


{Nj, i emiddles4,5,6 of triangle Tt } 



A function of form of the sub-families (203) are shown on 
figure ( 17 ) 



17.1 


17.2 


Figure 1 7 


It is obvious that N i 2 0 ifj, i e5 h (figure 17 . 2 ), but it may be 
observed that on figure 17 «* that « ,^*23 may take on negative val- 
ues. i Cc h 


In this case | we shall substitute for (198) the con- 
straints of inequality (204) ” " 


Qi 


n, 


J*?n! dxsO VN.€/»i 23 m! - Max (0,N.) 
n on 1 i n 1 1 


(204) 


Qi 


n, 


7( ^h^ > oh ) * 7N i dxS 0 V V*h 


456 


One way of satisfying them is to add to the criterion ( 193 ) the 
penalty functional ( 2 O 5 ) 


p . y I (jt ! 2 + 

‘ h 1 1-3 1 . 456 

i < *l - 3 u*” 


in 


+ .2 


(205) 


In the numerical tests, the second term of (:?Oj) shall practic- 
ally be sufficient. In tho case of approximation ( 102) po , the dis- 
crete penalty term of (20b) 




. f 

J n 


-•* * 1 * 
.u* n j 


& 


( 206 ) 


/6k 


A contribution of S due to 
a discontinuity of spued 
of positive intensity is 
shown on figure 18. 


U local 



Figure 1 9 


In the case k»2 the dis- 
crete discontinuity of 
figure 18 is calculated 
at intersection 4, 5 or 6 
in the direction of the_ 
bars of two triangles T** 
and t + c 7C. The illumina- 
« tion ♦ and — 
occurs by using the local 
speed on figure 19* 


1 .■* 

u 56 " “2 ' u 5 +u 6^ calculated in the middle of bar (23) view of 
and in the middle of bar ( 13 ) view of T, makes it possible to il- 
luminate Tj , T 2 in T" and T + in the following direction 

is exiting at node 5 =#> Tj -► T 

- u-, i® entering T 2 at node 6 T 9 -► T h 

In this semi-node, the discrete constraint to be satisfied is 
expressed by ( 207 ) 


56 


S i - Cu'nlj - (^0* - ^j)*nS0 (2 ° 1 * * * S * 7) 

discrete analogue of (206) may then be expressed in (208) 
^456^' * (Bi> 1 * lenffth of bar “ B i (208) 


if N. 


,-hich permits the final approximation of 102 in (209) to be given /6 5 


( 209 ) 


1" <W * V W 

♦h 


Note : it is possible to form a model when k=l . the condition 
of ontropy by adding a penalty to a functional odd power of a 

positive step of speed, which is impossible according to fluid dy- 
nnmics (decompression shocks) given in (210) 


l 


iCt-nf a j 


(be 


witha-3 


( 210 ) 


The discrete analogue is expressed then 


r +2 

S> 3 h “ 2. a R* HB.) with this time 

ie>v 1 1 

h 


(211 ) 


Numerical results using ( 2 1 1 ) with k=l shall be presented later 

on* 


M.O. BRISTEAU (20) may be consulted for the numerical approx- 
imation of the constraint of entropy by the artificial viscosity. 

9*2*1*4. Approximation of the Cost Function Gradient and of the Pen- 
alty PYmctional Gradient 


( 212 ) 


The 


cost function gradient <J'(4>),N.> ■ J| is approached by 


<J h ( V 'V • L pW h ,? 5 h ?,) i dx - H„ pw h )'' l/a<? V^i} <f V^i ,dx ( 212 ) 


v,: i‘ v h ' ?h <v h ' V v hg 


The discreto analogue of the menaltv functional gradient (202) 
is expressed by (i’ll) (21^) (21 r >) 


(213) 


<P i23*V 



2 . £. 123^1 with 

icN, J J 
a 

<V\ h >^v* )+ 



(i). given in (214)(2I5) 


6Qj(i) 


$ot*7N dx 
ft J i 


(214) 

(215) 


If the entropy constraint modelling (21 l) is used, we obtain /66 
the differentiation formulas ( 216 ) ( 217 ) 

<S 3h •»!> - 2 ,J 456 R* «Rj«)t( Bj ) (216) 

with fiR* given in ( 21 7 ) 

« R J(i>3{(?^ -V4>7)»n} 2 .{ (VN?-VN7)*n} (217) 

9*2.2. Case of Lifting Profiles (Airfoil Sections) 

9.2. 2,1. Approximation of Spaces V , V and 

g o 

If V and Vgh designate the approximations finite dimensions 
of spaces v and V g $ if V c is the sub-space of jj l (ft) defined in 7.3 
so that v ‘ ‘ 1 

V„ « (4> e H 1 (ft), <M trailing . $1 + - <f,J « £ , i any 
e • edge 'C T1 C value 


where C designates a cut in the domain occupied by the fluid exiting 
the trailing edge and joining a point of (Figure ?0) and that V* 
designates the approximation in finite dimension of space u 



{0 h U h cc°(n>, VT£t v f h l B p- 0 




i) 




wmmm 


Figure 20 


then the discrete analogue of (l 09 ) is given by (218) 

V ^Ph * **Rh e V h®^Ch 

9*2. 2.2. Approximation of the State Equation ( 1 1 2 ) 


In (218) ^Rh and ^NPh are the approached solutions ( 2 1 9 ) (220) 
of the variational equations (52) ( 1 1 2 ) 

I • dx * 0 V CO € V* • (ft £ w ^ /pi 

'& Rh h h V Ch ,<P Rh CV Ch (Z1 

The step condition on C is treated as a condition of pseudo- 
peri licity. We define £ which pproaches £ as solution of the 
discrete equation (220) h 

L VC h \ dX = f 0 P h ( ^h ) ^£h*H dX + [ p Vh dr 


V CJ. € V. 





9«2.2.3. Approximation of the Joukowskl Condition 


T «■ 

If Tpp and TgF designate respectively the last element at the 
extrados fresp, at the intrados) attached to the airfoil following a 
side of a triangle, and to the trailing edge shown on figure 21. 


3 



£ 


Figure 21 


we approach JK(£) by defined in (221 ) 

jyt) - I^h 1 ^- < 221 > 

If k=l . is constant on each trangle and the Joukowskl con- 

dition cannot be applied punctially on the airfoil, but only as an 
average on the two triangles T + - T”. 

If k=2 . •i* 2 linear on each triangle, we may selecte one of 

the nodes ( 1-4-2) oi (1-6-3) on the body or an interpolation of 
these points (refer to MARTIN ( 23 )). 

9.3. 3-D Flows 

The numerical implementation of the tridimensional flows is 
developed in detail in J. PERIAUX (9). In the case of lifting flows 
(for example, around a wing), a sheet of discontinuity (ND) must be 
introduced, originating at the line of the trailing edge (LBF) and 
joining r as on figure 22. 



P±S 





emplanture 
socket 


Figure 22 


* * 0 * quality 


The discrete Joukowski condition is applied on (LBF), trailing 
edge line of the wing* If Z(y) designates the circulation spread, 
starting from the socket of the wing up to the "pig" the discreti- 
zation of (LBF) into NBF points on which the sheet of discontinuity 
is applied, provides NBF Joukowski conditions (non linear) (222) 

nr ro \ _ I#* I 2 i A* I 2 


J W = iKi + - iKi - • k=i * nbf 


with decomposed as follows in ( 223 ) 

NBF 

*h * *NPh + l l ¥L ^Rh,K 
NBF l 

v v h®a *<£»,.„> 


. . v K defined 0241 
with v (nd)k in U24; 


* » ‘Wt 6 P k V Te ' t 'h 


■ i 

V 





If Tt is a tetrahedron of ft then ND shall be the trace of /69 

tetrahedrons having at least one node belonging to the sheet of 
discontinuity 

(ND) + shall bo the trace of the sheet , view from above 
(ND)” shall be the trace of the sheet, view from below 


The definition + and - are defined from the line of the trail- 
ing edge by designating by + the extrados of the wing and by - the 
intrados of the wing. I shall be observed, then, that the discreti- 
zation of (ND) h is composed of a set of triangles (see figure 22), 


O.C, ZIENIEWICS (2k) may be consulted for the approximation 
k=l, 2 and the coordinates of surface area (L^) used in (225) as 
well as the derivatives 



4 ^ 10 

i, v-i * v wv ! 5 i T c? k 


(225) 


of the functions of forms appearing in the exact integrations. 

10. - MIXED APPROXIMATION BY THE METHOD OF CONFORM FINITE ELEMENTS 
OF THE NAVI ER- STOKES EQUATIONS 

10.1. Summary 

This chapter presents the mixed approximation of the Stokes 
equations and the Navier-Stokes equations by the method of conform 
finite elements taken into consideration in chapter 8 . For simpli- 
city ft shall be assumed to be a bound polygonal of R , but the nu- 
merical implementation extends to the domains of r 3, the applica- 
tions of which shall be presented during the presentation of numer- 
ical results in chapter 12 , 

10.2. Appr ox imation of the Functional Spaces 

If ^ designates a standard triangulation of the domain & , the 
following spaces of finite dimension ( 226 ) ( 227 ) (228) ( 229 ) ( 230 ) 

( 231 ) shall be used subsequently /70 


1 


H 


oh 



| T € 


" “J - <v 


(v. e 0°(ST)) 2 , t 


h 1 T 


l v v h • *hW 


1 

, * h l T - 0} 

« <p,) 2 , v ic tr h ) 


(226) 

(227) 

(228) 


ss 


(229) 


with 2^, an appropriate approximation of z. 

"sh - : < vV ' \h 11 “oh • j a \-\ “* • /Hv* v 


VV 


A widely used variant in numerical tests consists of 
defined in (231) 

\/2 ■ <V2* (<J °® >2 - V2ll £ ‘ P l> 2 - V T ‘ & h/2> 


(230) 


(231) 


where *V2 is the triangulation obtained from bysubdivision of 
each triangle Tc into 4 sub-triangles obtained on figure 23 by 
joining the middles of the sides* 




Figure 23 


10 , 3 * Approximation of the Steady Navier-Stokes Equations 

The approximation of equations ( 1 32 ) by the mixed method ( 1 87 ) 
is given in (232) 


111 


<v 


Find W zh 80 that 

j n dx + j ( v^v ( v\ )dx ■ 0 » v { \> k u w oh (232) 

We shall find in P. LE TALLEC (26) reasonable assumptions on 

and v so that (I’h) permits one solution. Moreover, passing to 

the boundary is the solution of problem ( 1 32 ) 

^0 0 { “» o} 


(233) 


the 


It may be observed that if *^h " ® is set in and in (P^), 

Taylor-Hood (25) scheme is recovered 

10. *l. Mixed Approximation of the Unsteady Navier-Stokes Equations 

Subsequently, k ■ At shall stand for the discretization time 
step. Presented now are two possible discretization schemes of 
which one is semi-implicit and the other one is entirely implicit. 


10,4.1. Semi-implicit Scheme 

This is the discrete version of scheme ( 1 52 ) ( 1 5 3 ) • It is com- 
posed of (234) (235) 


-►o 


u^c , is an approximation of u ° f given 


(234) 


Then, for 
r-^n+I ,n+l-. 
{u h -*h J 


n 2 0, by using ( 232 ), we obtain (235) 
from -* n by solving ( 236 ) 

“h 


**-n+ 1 _-*n 

l n “nr-^h dx * - -f n <\^>%-<v\ >dx 


(236) 


v l »h-V tW oh * v 


It may be noted that ( 236 ) is a sequence of discrete Stokes 
pseudo-problems and that the scheme ( 235 ) ( 236 ) is a truncation 
error 0(At) a nd on ly conditionally stable . 

10.4.2. Implicit Scheme 

The scheme taken into consideration above in (237) (238) is an ZI2 
entirely implicit two step Crank-Nich olson scheme 

•>0 -*-1 

u >Uh given ( 237 ) 


Then for n ~ we obtain by using (232) u? + * from * 

the solution of ( 238 ) “ 


3 -Ki+1 _->-n 1 -t-n-1 

2 % 2u h + 2 % 




•v, dx + v 
n 


dx * f (tr'.fts?* 1 . 
ti h Jn h 


( Vv« h )dx ■ 0 v 


(238) 


It may be noted that ( 238 ) is a sequence of discrete non- 
linear problems analogous to ( 232 ) and that the scheme (237) (238) 
has a truncation error o(At^) and uncon< *itionallv stable . 

10,5* Least Squares Solution of Discrete Unst eady Navier-Stokes 
Equations 

10.5.1, Discrete Mixed Formulation of the Problem pj 

We are taking into consideration in this paragraph the discrete 
analogue of chapter 8,2.3. » the mixed formulation of which is a gen- 
eralization (239) of ( 232 ) 


Find € W zh 80 that 

otj u^'v^dx + vj " 

■ L f oh‘V x * { Q f ih -( v \ >d * v( Vv £ w oh 


(239) (Pp 


Two terms may be observed in (239) 

4, corresponding to the choices of the quantification 
“ * ^ time sbheme 


k 


(2 


(2uJJ - j^ _1 )dx ( in (238)) 


- f t . density of external forces 

I n 

Ph is a nonlinear problem. 

10.5*2. Least Squares Method of P^, 

By analogy with ( 1 58 ) of chapter 8.2.3* » the least squares me- 
thod of Pj!j given in (240) (24l) (242) is taken into consideration 


in 


fv u J h (V h’V vlth 

tv h*V eW zh 


(240) 


VW • I f n I V^hl 2<lx + 7 | n ^ ( W* 2{lx 


(241) 


where t is a function of {v, ,<{>.} via the discrete state equation 
(242) “ h h 

“la Vv* * v / £1 ^h*^h dx - J n ? 0 K T ’h dx * f/i£<vV dlt 

■j'V^VvK 1 '* ■ v ( V“h )£ “oh * <W £ “*h 


( 242 ) 


10 »5.3o Calculation of Gradient J. 1 

n 


The differentiation of criterion (24l) is given in ( 243 ) 

5J h " °/ fl ( v^h ),6( W dx + v J * ( VV‘™ ( vV dx 

V <«V 6 V * W oh 


( 243 ) 


whereas the one of the state equation is given in (244) 

i { v\ dx * 

- I <V ?)d V<V ? “h )d * 


(244) 


uith 

Since ^h^h’^h^h^ W 0 h’ it: is P ossible to express the variation 

of criterion gj uniquely as a function of £ by using (244). 

“ h 

We obtain (245) by selecting {q * ^ jf 

h h ’h ’ “h e ^h' x h } 

i ■ ./ «.*», . . /s 

* ?M h -x h ))dx - / !j <v h ^)^ h .«7 h -| h ,4 Mh - Xh))djI (245) 

By putting (243) (245) together, SJ h is finally given by (246) {jk 

6J h “ °L ( V^h ),(S \ dX + V [ ^(v h -1u)*^V. dx 

J n h J« hh h (246) 

+ J« 6 V ? >\ * ‘V^^vV^h-Xh^* 

By expressing that 5 ^ -^<J£(v h .$ h ) , (n h ,u> h }> , ma y be identi- 


fied with the linear form ;W Qh ^ R defined by ( 247 ) 

<J h ( vV- t \-“h )> ■ a f<v^\ dx * l ? <v*h>-^h d * 


'fl 


a 


.n'W * *V*)\»-«vt h )*?« h -x b ))^ 


(2^7) 




r 


10.5.4. Conjugate Gradient Solution of (240) (24l) (242) 

The algorithm given above ia the discrete analogue of the one 
described in (l6o),.,(l65) in chapter 8.2.4. It consits of 3 
phases. 

Phase^O t Initialization (248) (249) (250) 
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( v<> € W oh 8 lven 


(248) 


Calculate (8^,0?) solution of the discrete variational equation 

(249) h h 


a L 3-v»* ♦ 


}> 


v {i V“h )£W oh • { ^- e h )e “ 


oh 


(249) 


with defined in (247) 


.et 


<V 6 h°> 


(250) 


Phase 1 : Descent (251 ) (252) (253) 


A" - arg Jin 

•*m + ! -*-m , m*m 

U h " u h ' X % 


,m+l 

1*1 


Cv - iV 

n h 


(251) 

(252) 

(253) 


Ebase_2 t Construction of the New Direction of Descent 

Compute {t n+ l Q nj+| i solution of the discrete variational equa- 
tion (254) 8 h » 6 h } 

i sr'v * v Lk‘-k <w> 


* { vv< «oh • w oh 

Then calculate the coefficient of conjugation y®** in (255) 


(254) 


_ °/ n ■> 

°/fl djt + v [ l^s“l 2 d* 


(255) 


The new direction of descent is given, then, in ( 256 ) (257) 


hj* 1 - IJ ♦ 


(256) 


6® + Y m+! T m 


(257) 


Do m=m+l and go in (25l) 

It may be observed that each iteration of the algorithm (248)- 
(257) requires the solution of several discrete Stokes problems 

-one Stokes problem to solve the state (242) {^ m+ l ^^m+1 j 


with 




-one Stokes problem to calculate ( gjf* * , G®* * } from {^® + * 1 } 

, n h n ,r h 

and /cm+1 v m+|, . ..... 

is>h »X h i via (254) 

-several Stokes problems (= 3) to calculate 

Flow chart 1 of the unsteady Navier-Stokes algorithm is presen- 
ted below, The discrete solution of problems g is presented in 
the follow. : chapter, 10,6, 

10,6, The Stokes Algorithm (Discrete Case) 

10,6,1, Introduction 

The presentation of the Navier-Stokes equations (discrete in 
tho steady case (it*. 3 ) and unsteady caso (l0,4) in the form of a re- 
petitive sequence of discrete Stokes problems, implies a highly ef- 
ficient numerical algorithm of problem ) ( 258 ), 1 “* 


SOLUTION 2-0 BY THE CONJUGATE GRADIENT OF 
THE LEAST SQUARES METHOD 

0 Initialization : selectu* 6 W£ 

* COMPUTE < g *• , 1) > a ( ) V ~t) € W, 

r.ewo 


1 Cl leu 1 atoll , solution of the 

k State 



3 Dascint : Catcul* X n 


^ n e Arg Min j(v!) + A)T") 
% > 0 * * 

set v , n + A" If " 

k k k 


21 Oirichltti 



A Construction of the' nev direction of descent 


Set ]f" +, « 1 %*'+ hj 


k '.index (time loop) 
n index (control loop) 


.7 C icHliti 


3 (fichlits 
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* Find ( * ^zh 80 ***at 


<«•><**) 


/ Q VV^f/vK ix ■ l a \A ix * | fl i ,h<V*h> <lx * <Vh> « w oh ■ 


til 


In the following presentation, we shall find again the discrete 

analogue of 8,3 in the solution of ($ ) 

ah 

( 27 ) (28) may be referred to for demonstrations of the theorems 

used. 

10.6.2. Characterisation of the Solution 

We may easily prove that ( 238 ) has a unique solution which is 
the one of the problem of minimization ( 259 ) with distributed linear 
constraints 


, Min Cf 

t(v h'*h )><W zh 


,+• I 2, 

I V dx ♦ 
n 




n 


£ *v.dx 
oh h 


- [ 0 W* 


h )dx (259) 


where it is recalled that 

M .h • { 'VV‘ v zh' < “ih' f/v ? 


V x ■ 


n 


* -v h v* v v “i» 


The number of constr ints of (259) is dim(Hj,). We may combine 
with (259) the Lagrangien ^ • V * h! » h 1 * ® defined by (260) 

h " h no n 


■ VvV * | jv K dx '{/' Vh dx 


( 260 ) 


where ih^h’^h' equals ( 261 ) 


VW ■ 2 




|v h l 2 nx. 


if l» 

1 h. 


2 f - - r 


fv. I dx - f .*v.dx 
h' Jfl oh h 


Q 


? ih ( v\ )dx ^ 261 ^ 


( 259 ) being a problem of minimization with linear constraints of 
finite dimension for which there is a solution and a distributed 


Lagrange multiplier 


>h €H h 


so that 


... 1 is a saddle point of 


h on V 2 h * K oh XH h vith { W solution of (258) (259) * 

The extronio conditions of ^ at point t , p^} (262) (263) 

(26'i) characterize the solution of (258) 


( 262 ) 


I a ****** ' |/ih-^h 1111 v V < h . p„« < 

L Vh- * */«. V^VV-jA**^ 


dx 


V. -r V . ■* ,, 

h ' °h ’ “h ^ V 2 h 


L ** - (k < h <« * v < 


(263) 


(264) 


From (262), it may be deduced that the Lagrange P h multiplier 

is the discrete pressure. 
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10.6,3, The Space ^ 


By using the observations made in 10,2,, l h is introduced as a 
supplement to in , i,e, 

"J • H oh®*h 


with ■ dim Cfyj) • 

In practice, by using the Lagrange finite elements, //?. shall be 
defined as follows ( 265 ) ** 

V\ y h ! T - 0 V Tc X2 h so thr.t3Tn 3Q - 0 (l. 


*h has a finite dimension N^, It is the number of nodes of 
belonging to ^ Moreover, ( 265 ) implies that “ 

supp(y. ) T7-, ■ l_J T with lim mes (HI. ) ■ 0. 

h 1 h Tn3(W h*0 h 

is shown on figure 24, 
h 

10,6,4, Converting Problem J fth into a Variational Problem ( ) in 
10,6,4,1, Approximation of a(*,*) 


Tn roforoiico to tho ebsorvations made in paragraph S.3,'3,, tliat 
if V is sufficiently steady, the Green formula loads to 



f 3lK r . 

a(X,u) ’ 'j r *T V iT " j n dx - j w x 

= Jn dx + f y dx = -f (fy 

n 'ft A Jn 


A +U A^ dx 


where m is a measurement of P in £ 2 . 

Now let ■ 5 l h f ' J h G ^h' If we define a h^* ’ ^ 
sequence of problems (267) (268) (269) {270) 


• ?p Xh -?q h dx = 0 v q h £ H oh • p \h'V H oh 

4 4b 

f -*■ ♦ . f ^ ^ ( ~ ■*" 

D II. »V Ujf + .U._ • '.'V. ux = - 


u. 


v n 


v, dx V v, € V . , u,, e V . 


f 2f 











* {Mf* i , t 1 - 

,. _• . -. - j .,' .'" 



Then, the theorem (10.6.4.1.) demonstrated in (30* dis- 

crete analogue of the theorem 8.3>3»1«» characterizes the properties 
of the bilinear form a h (.,.), 

v t« tr 

Theorem 10.6.4,1, j Let us assume that "h* T has at the most 

one side c 3ft , therefore ajj(.,,) is a bilinear, symmetrical form 
and is defined positive on 

C^/Rh) * where 

**h * 1fJ h = cste on 3ft} 

Based on theorem 10.6.4. 1.^, we can now convert the problem S ah 
into a variational problem in ( h thanks to theorem 10.6.4.2., 
discrete analogue of theorem 8.3.3«2. 


10.6.4.2, Approximation of (e) 

A 

Theorem 10.6.4.2. j Let P. be the discrete pressure and h the trace 
of on Therefore if theorem 10.6.4.1. is verified, then A^ 

is the unique solution in ** /r of the variational linear problem 
(i^) (271) h h 


X h e V*h 


“h^h^h* = 


8 


^oh- oh )-^ h 


(271) (y 


P h 1 U , Ip 

where ofl °h are defined by the sequence o^ nroblems (2J2) 

(273) (274) 


80 


O ?!> oh - K dx ■ L/'»' fq h d * 


V V H ih- »oh e H oh 


(272) 










. 1 .jiij m rv! j l.., i j r : vr^r.T 






(273) 

(274) 


10,6.5. The Solution of Problem (E^) 

10.6.5.1. Summary 

The choice of the method used for solving the problem depends 
uniquely on industrial applications. For 2-D fluid flows, the num- 
ber of boundary points N h (*100) with dim N’ h «dim H^, the solution 

of (E ) by a direct method is preferred, for the core space and man- 
ufacturing time required for matrice is relatively compatible 
with the current size of large computers (370/168). On the other 
hand, for three dimensional applications ( separated , flows around a 
wing with high incidence), the number of boundary points N. (- 1000) 
results in an unallowable core use and computation time and in this 
case, a conjugate gradient type iterative method is preferred for 
the solution of (E^), which does not require information about 

Both methods are expanded in the following text. 

10.6.5.2. Solution of (e^) by a Direct Method 

10.6.5.2.1. Construction of a Linear Syst em Eq uivalent to (Ejj) 



General s 


tjy 

The space /( h defined in (26'j) being of finite dimension, let 




, a base of 77?^* That means that y y e 77!^ 


N. 


1 * N 

y h = J, Vi 5 Vh = { V-‘-\ e]R h ) 


(275) 


The functions w i are defined as follows : 
V i=l .... ,N h 

V V h w i CP i) = 1 

w £(Q) 3 0 ¥ Q node of*^ , Q t P. 

a •> 


/B 1 


( 276 ) 


I 

j 


J 

•i 






J 


'll ■ Tin' ifl’ItfH 1 " ii i 


Pi 


///f. JupparT 

Da 



Figure 25 


The hachure zone of figure 25 
represents the support of w^. 

With the definition (276) of the 
w if that means that in ( 275 ) (P^) . 

The problem ( E. ) is therefore equiva- 
lent to the linear system ( 277 ) 


*h 


- J n (^ oh +u oh )-^. dx , lsisN h . 

Set Yj * WV J \ - < a £j>i<i,j < :1 ; ^ - \ ^ 0 h + “oh^ w i dx ( 2 77) 


b. - 

n i i«| 


According to theorem 10.6.4.1 . , the matrice A. is complete, sy- 
mmetrical and semi-defined positive “ 

Construction c f A : A^ is constructed column by column according to 

the relationship a^ = ah(wj» v i).To compute the column of A h , 

the sequence of problems ( 267 ) • o . ( 269 ) is solved for \ = w j and 
( a i j )i=l, . . .N h is deduced by using ( 270 ). Each column of A re- 

1 h / 

quires, then, the solution of 4 discrete Dirichlet problems (5 in 
the case of ftc • As the matrice A^ is symmetrical, the problem 
may be limited to indices i > j . 

Taking into account the choice of ‘h , the integrals defining 
( 270 ) involve only the functions having a support of about ^ 

(Figure 25 ). 0 '“ 

Flow Chart 2 of the construction of operator A is presented be- 

lOWo 

Construction of b^ s to construct the second member of ( 277 ), the 

sequence of problems ( 272 ) ( 273 ) (274) is solved, which requires 4 

discrete Dirichlet problems , = . __ ,,3„ 

( 5 1 f ft c R ) . 

Considering the choice of^h, the integrals defining the sec- 
ond member of ( 277 ) involve only the functions having a support in 
the proximity of (Figure 25 ) 0 


10.6. 5. 2, 2. Solution of System A^A^ = 5 ^ by the CI 1 I 0 ski Method 




A) COHSTP UCTiOH OE t'UPEHATEUQ SyMETQIQUE 
OEFidl POSITIF A (i, j ) CO LON N £ PAR COIQNNE 

KEY A) CONSTRUCTION CF THE SYMMETRICAL OPERATOR DEFINED 
POSITIVE AT (i,j) COLUMN BY COLUMN 
B) For all nodes... 

s)r — ► Pour tou« lu noeudi Ief 


.V"m dn = 0 
Pk If = 5 IK 


^«(u'l ( .fr+v^u k .^Ndii = / v’ p k .lT d 

u k |r = 0 


Jsi m dA = k M d,a 

4*k lr = 0 


a(i,o)= k.7 k Mj + 

d A 

support J , J 6 T 
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Account taken of theorem 10.6.4. 1. and of the definition of 


■&. 


Ker(A h ) - {veR^Vj ■ v 2 " V N.^ 


and since the matrix Ah is singular, it is necessary to fix a com- 
ponent of A^ ■ 0 , for example) in order for the sub matrix 
h 

a . (a..)i • ;<n -1 to be symmetrical, defined positive , 
r, 1 j 1 1 » j- 

The sub-system to be solved is therefore expressed ( 278 ) 


Vh x h • b h uh,r * 


( 278 ) 


r h X h * * X 1 * b h " {b l ’ b 2’ ,,,,b N, 

n n 


( 278 ) is solved by the Choleski method via the standard factori- 
zation ( 279 ) 

* ^h b h where b h Ls a non singular lower triangu- (279) 
lar matrix 


In summary, the solutions to be computed to obtain the solu- 
tion / 1. - \ derived from (E, ) by the Choleski method are the fol- 

lowingVV h 

-4 discrete Dirichlot problems to calculate P oh ,U oh^oh b h 
(5 if n c R 3 ) 

-4(N -l) discrete Dirichlet problems to construct \ 

^ (5(N h -l) if Q c R 3 ) 

- 2 triangular or descent-climb systems to calculate A^ 

HA " b h ! Ltf h x h ■ 5 h 

- 3 discrete Dirichlet problems to obtain p^ and from A^ 

(4 if c R 3 ) . 

Flow Chart 3 of the rapid Stokes algorithm is presented below* 

In practice, the matrices of the Dirichlet problems are fact- 
orized once and for all outside of the control loop. They are two 
symmetrical matrices defined positive, one approaching -a by ele- 
ments Pi, the other ald-A by elements P2 ( or Pf on a triangulation 
-f? , defined in ( 213 )» 

ti/2 

10.6.5.3. Solution of (Ej 1 ) by a Conjugate Gradient Methoc 
Gerxoral : It 3 - s interesting to solve (E ) by an iterative method 
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rapid stokes algorithm 
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which does n t require the EXPLICIT computation of Ajj, but only, at 

each iteration, the solution of 4 discrete Dirichlet problems 
(5 ifftcR^), It is subsequently interesting to introduce the iso- 
morphism 

N h 

r h : ^ defined by 


Vh ' ( V U 2 *(,> *V\ 

"h h 

with^h " | ^i w i ’ w i e ^h dase ^ ncrodu CQd - * n ( 27 -i) 


(•f«)h sha11 designate the standard euclidian scalar product in 
N the corresponding standard, 

R h * IMI h 

By using the observations above, the two members of problem 
(Ejj) ( 271 ) are expressed 


WV ■ ^hVVVh v X h ,u h c 

f a dx " ‘WiA ¥ V *>h 


( 280 ) 


Description of the Algorithm in (281 ) ( 82) ( 283 ) (284) 


EbaSfl-Q-A Initialization 


r l° tA 

hh £ 1S given arbitrarily 

8 h = A h r h A h~ b h 
.o o 
h h - 


(281 ) 


Then, for n>0, being known, compute^ ' ,g“' 1 .h^' by 


A. n+, . & n+I un+l 


Phase 1 : Descent 


(h n ,g n ) ( 


IM. 


<v>;>h 

Vn - p h ^ 




(282) 


ji+1 n 


8 h *«h-Wn 


( 283 ) 




£hase^2 t Construction of the New Direction of Descent 

II „ n+ I ll 2 


" I* e hHh 

C ■ C » Y h h„" 


( 284 ) 


n*n+l, go in 282. 


Notes t 


As the matrix A, is symmetrical , semi-defined positve, it may- 
be shown that the sequence {X^} converges toward ^solution of 
(E^) c The component of X^ defining the pressure level is the same 


one as the initial pressure 


y 


The implementation of (28l)..« 


(284) requires^the solution of 4 Discrete Dirichlet problems to 


obtain p 
to compute 


h.h n ’ u n.h n ’ ^n,h n 


at each iteration (5 if Q c it) in order 


A, h via (285) n ,, . , A .n _ x 

" ( a h = (A h h n* h y h ) h 


I 


( \,h" * “h.h^’^h d * 


(285) 




The pre factor! nation phase of the discrete Dirichlet matrices y 
recommended in the direct method, is also obvious, upstream of al- 
gorithm ( 281 ) . • • (284) and leading to considerable gain in calcula- 
tion time. 

10. 6 . 5*4. Acceleration of Algorithm ( 281 ) . . . ( 284) by Preconditioning 

Let h 'h "‘h be a symmetrical bilinear form defined positive 
to which the symmetrical matrix defined positive S h is related via 
(286) 


a h (X h’V 


(S h r h X h 


•Wh 


(286) 


S. is an auxiliary preconditioning operator in the sense of 0. /87 

AXELSS0N (32). The conjugate gradient variant using a scalar pro- 


duct 

( 288 ) 


(( VWs 

(289) (290). 


< x h* s h y 


relating to S h is defined by ( 287 ) 


0 * 


Chase 0 : Initialization 


(287) 


N. 

r.XfcK selected arbitrarily 
n n 

o 
o 


o A ,o 

«h ■ Vh x h • b h 


c“I O 
h 2 * * * * * 8 * h 

For "iO X”, g”, h£ kD'--«, compute l£ +l , g” +l , h£ +1 by 
Phase_l : Descent 

/ l. n n . 

<VhX> 


r X n+I 

r h A h 

n+1 


i - p \ h h 


( 288 ) 


Phase_2 > Construction of the New Direction of Descent ( 289 ) 

, n+1 C -1 n+1. 
v . (S h ,S h 8 h 

0 fr , 11 C _1 

8 h ,S h 8 h^h 

u n+l _ c -!„n+l , n 

h h S h 8 h + Y n h h 

n*n+l, go to ( 288 ) 


(289) 


Notes : If S. a Id (identity matrix) is selected, algorithm (281... 
(28 ^ ) is found again. 


Different choices of S h are proposed by GLOWINSKI-PIRONNEAU (29 
( 29 ) guided by two c .fferent types of contradictory arguments (info- 
rmatics and theoretical). 


1. Select S. (.,.) leading to a hollow or even diagonal matrix 


In 


Shi 


s case, 


S, may be factorizd once and for all by 
h 


the Choleski method S = T upstream of the algorithm 
/ fa h h 

(informatics argument). 


2. Since a h («»«) is an approximation of a(,,.) defined on 

H _l/2 (r) and elliptical h"’ /2 (r) » select S ( • » • ) approxima 

tion of a bilinear form S(.,.) also elliptical h 1 ^ “ (r ) 

This alternative, however, leads to a complete matrix S 

(theoretical argument) 

We give to ( 290 ) (291 ) (29 ) three possible S},(,,,) leading to 

sparse matrices, provided that the boundary nodes have been 

numbered properly (minimum band width). 





I 


vw 


s h ( W 


j r Vh « 

/ n x h»h ® 

K-K * ■ 


(290) 

(291) 

(292) 


f 


Assuming defined in (275) ( 276 ) and that the Lagrange finite 
elements are used for the problem ( 271 ), it is then possible through 

numerical integration to combine with ( 290 ) ( 291 ^ .ilinear forms for 
which S h is diagonal* This is the approximation f 293) for (290) 


<Vi. 


IA 


5 h ( W 



described on figure 26 


(293) 



Whereas approximation (294 
294) is related to ( 29 1 ) ? 
describes figure 27 


S h^h»V " | 3 mes (supp (M.))A.y. 


'////A Support 



Mi 


(294) 


Figure 27 


Hgfi 




In (294) we recognize the scalar product L approached (,,.) 


defined on V h by ( 295 ) 


N+l 


^ f h ,g h^h"(N+l)_ I f h (M i ) 8h (M i ) if£lc rN ’ N * 2 or3 

T e i-| 1 " 


(295) 


with mes 


l? i 


area or volume of T 

nodes of triangle or tetrahedron T. 


We shall check whether the matrix S i 8 diagonal by using the 
definition of B h give.i in (275) ( 276 ) alid (294), In fact, in this 


case 


r, hWh • J oe«(*upp MpXj. 

Finally, it seems interesting to select S. the inverse of 
the matrix (292) in an approached space of ^-1/2^" f i.e. ( 296 ) 


-I 


(X h*V " dx 


( 296 ) 


Various numerical t9sts of possible conditioning of (290) (291 ) 
( 292 ) have been applied to the solution of the discrete Stokes pro- 
blem via ( 287 ) ( 288 ) (289) e Tho rapidity of convergence (number of 
iterations) and the calcu; ation time are presented in chapter 12, 

11. - ON THE METHODS OF INCOMPLETE FACTORIZATION 

11.1. Summary 

This chapter deals with the difficulties of informatics imple- 
mentation of least squares algorithms on two and three dimensional 
configurations of large dimension . 

We show how to use the methods of incomplete factorization as 
auxiliary operators of preconditioning or as auxiliary metrics in 
order to overcome excessive transfers of data on auxiliary memories 
(disk and or bands) outside of the main computer center. 

11.2. Auxiliary Operator of a Problem of Model 'volution 

We shall now consider the parabolic linear problem of standard 
evolution defines in ( 297 ) ( 298 ) ( 299 ) 



f(x,t) in fix ]0,T[ 


(297) 


0 on f x ] 0,T [ 


( 298 ) 


$(x,0) ■ $ Q (x) when t»0 


(299) 


Z22 

where tJ designates a bound domain of R n of boundary p with f and 
$ 0 sufficiently stable. 

Any quantification in implicit time of (297) such that 


^ $ k+I - A$ k+I - ~ $ k ♦ f(£,kAt) ( 300 ) 

|( •+ s 

with $ - 0(x,kAt), At time step leading to the solution of a lin- 
ear system (301) after quantification of space (of finite differ- 
ences or finite elements type) 


A$ 


k+l 



(301 ) 


A* (d .) lii i^N 

where ij is usa ny a positive defined symmetrical ma- 

trix (N X N) with half band width m (N representing the number of 
nodes strictly included in the quantified domain . 

Since (301 ) must be solved numerous times and that A is inde- 
pendent from k t it is better to use a Choleski type direct method. 
Since A is symmetrical, defined positive, there is an inversxble 
and unique lower triangular matrix L, having the same band width 
as in , so that 


A « LL fc 


(302) 


with > 0 I 1 £ i < N 

ff . 

where ii » 1 are elements of the diagonal of L. 

if jj, are elements of L so that 

ij 

*£ j ■ 0 if 1 < i < j s N 


We bring back the algorithm of factorization of A ( 303 ) ( 30** ) 
f for j*l 




'I 1 


l l 1 


11 ; -n 


, V 2 < i < N 


( 303) 


(304.1) 


I 


Tor 2 S j < N 


j’ 1 


Ai . * (a. . - 7 j i.) 

jj JJ 4 


2J/2 


k-I 


jk' 




i-i 


ij ' *77 U ‘J ' J, V j* 1 S 4 £ K 


(304.2) 


Once L is calculated, the determination of ^k+l is immediate 
via a "descent-climb" ( 305 ) 


L 4> 


L t *ke| 


(305) 


In industrial applications N may be very large ! 10000), 

making the storage of A and L in the main core of the computer even 
impossible. Moreover, even though the non zero elements of A are 
not numerous (A is a sparse matrix); as i ^ the matrix L, it is un- 
fortunately always full . 

Consequently, auxiliary core stations (disks or bands) must 
therefore be used, and this requires costly data transfers, which 
becomes excessive in an industrial context (problems of input-out- 
put, process time, etc,,,). 


In order to preserve the advantages of direct methods such as 
the Choleski factorization, it is dosireable to find a sparse lower 
triangular matrix ^ close to L regarding their spectrum, and kept 
COMPLETELY in the main memory. 

With £ it is possible to construct A for (306) 


We substitute, then, 
7 .k+I 

A v 


** f. 

A « LL f 

for ( 301 ) the 
- (A-A) i k + F k 


iterative 


process 


{ 306 ) 


(307) 


(307) 


In ( 307 ) A plays tho role of auxiliary operator of A. It may 
be pointed out that the strategy to be adopted is different in sel- 
ecting ~ depending on whether ('S 07 ) must be solved once or :.v v era l 

times . ^ 




In the first case, we shall look for incomplete, fast and ef- / 92 

ficient factorizations in storage usage L (see MEIJERINK and VAN DER 
VORST (30)) or similar iterative techniques (see VARGA ( 31 ))* AXEL- 
SSON (32), MANTEUFEL ( 33 ) ) » whereas in the second case, it is worth- 
while to perform a significant computation upstream of (307) (proces 
time, memory) to benefit from extremely fast solutions at each At. 

■H 

According to AXELSSON (32), A may be used in another way, by 
having it play the role of a preconditioning matrix for a conjugate 
gradient solution of (301). If a is used to define the scalar pro- 
duct (308) in R N instead of the usual scalar product (309) 


«M> ■ A 

(308) 

($,l i>) * 

(309) 


Therefore, the conjugate gradient solution for solving (301) 
corresponding to the minimization (310) 

J($) 1 A 4 - F$ (310) 

is given in (31l) ( 312 ) ( 31 3> 


Eb&SS.Q * Let <j>° be selected arbitrarily 

Calculate g° = A$° - F , , 

v * (311) 

R° = A _, G° 

Se - H° - R° 

then, for n > 0 , assuming .n n „n as known, calculate <f, n+ l G n+ * 
H n+ 1 by ’ ’ ’ 


Phase 1 s descent 


t 

A n = Arg min J($ n -AH n ) = J iT ^ 
\>0 H n A H n 


0 n+i = $ n - AV. 


(312) 


Phase 2 t Construction of the New Direction of Descent 




G”* 1 - g” - a” ah” 
R"*' . A -1 G"* 1 

Ml G"* 1 ^”* 1 
Y ■ — 

G nt R n 

H n+I - R n+I + Y n+1 H n 


do n=n+l and go to ( 312 ) 

Note : The closer a is to A, fever the iterations are required 
to obtain the convergence of (311 ) (312) ( 31 3) • At the extreme, if 
= A, the algorithm converges easily in on iteration. The number 
of iterations required for convergence is a verification, a poster- 
iori . of the efficiency of a. "* 

11. 3« Auxiliary Metric Related to a Functional Least Squares Method 


A situation similar to 11.2 exists for another class of equa- 
tions with nonlinear partial derivatives : this is for solving tran 
sonic and Navier-Stokes equations expanded below by the functional 
least s uares method. 

We combine with (31*0 (315) 

tT(4>) - tf«p(|fy>| 2 )fy = f («) (3i*» 

<M r ■ 0 < r > (315 

where p is a nonlinear, positive, bound. . given value of jV<pj^ 

The minimization (316) in of (31**) 


win J( W = IICW) -f || * 

H ‘(ft) 

is equivalent to the optimal control problem ( 317 ) 

min { [ |Ve| 2 dx | Ae = C(<fr)“f ,‘e| r =* 0} 
(freH'ttl) J Q . i r u; 

o 


(316 


(317 


The qua * fication of (317) leads to the problem of minimiza- 
tion in wi constraints (318) 




min {E C BE | BE - T($) - F} 
$e* N 


(318) 
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where B designates the matrix corresponding to the discrete Dirich- 
let operator and T the transonic operator obtained by quantifies 
tion of ( 3 1 ^ ) • 

Now, let us assume we know how to construct B close to B as in 
11.2, then in place of (318) we propose to solve ( 319 ) 


min {E* B E|B E = T($) - F} 


(319) 


If ® is defined positive (318) and (319) are strictly equiva- 
lent. . However, if g is not selected well (319) may be not as well 
conditioned as (318) and consequently, a conjugate gradient solu- 
tion of ( 319 ) shall require considerably more iterations than of 
(318). 

In this case, ® defined in (320) is the auxiliary metric of the 
nonlinear operator T 




(320) 


11.4. Construction of the Auxiliary Operator A (resp. metric B) 


We shall expand, in this paragraph, a methodology giving access 
to a class of sparse matrices ^ or g close to A or B. 


Let A = (aij)lsi,j<N a positive defined symmetrical matrix 
with half band width m so that (Figure 28) 

“ « I- -| . (321) 

a.j = 0 if I i-j| > m 


Since A is factorized by the Choleski method ( 304 ) (305) 

A = LL C 


(322) 


L is a lower triangular matrix, also of band width m. Further- 
more, it may be observed that even if A has MANY zero elements IN- 
SIDE the band (Figure 28), it is not the case of L, which has NONE 
(Figure 29 ). 


Definition : Let us define in (323) the set of indices K of 
zero elements of A inside band m 


K = Ui.i) i 



(323) 


a* 




CHOI ESKY FACTOOIZATIOH 


^i»*** ( * 


A * ll' 
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Figure 28 
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and lot us designate by n(K) the number of elements K, Since a 
positive constant C is now given, it is possible to define 2 aux- 
iliary operators !_ and 7, as follows 

C L C 



is defined (324) 
by 

“ 0 if (i.j) e K & |*„| s C 

■ otherwise 


L£ is defined (325) 

f = 0 if (i,j)e K& |£ J <C min (£..,£..) 
) , J ^ i.j 11 JJ 

' ij * ^ij °'* : l 1Grw ise 


(324) 


(325) 


The constructions of 7 and l* bring to light the following 
observations • C C 


1 . 
2 . 

3.. 

4. 


If C ^ min |£..| then L_ = L 

’ J |A..| 

IfC£ " in 1 then L i ■ L 

1 9 * • 11 11 

1»J 

If C S max |S...! then ((i,j)|£.. = 0} = K 
i,j |£..| 

T f c ^ max } theil {(i,j)| V. 

n JJ J 

A > 


= 0} = K. 


In cases 3 and 4, 4,1.^ and L' have their non zero elements lo- 
cated in the same position as those belonging to A and are very 
close to the incomplete Choleski operators proposed by MEIJERINK- 
VAN DER VORST ( 30 ) and D. KERSHAW (34). Nevertheless, they con- 
struct L DURING the factorization of A (which means that L is NEVER 
constructed !) and economize store usage with the possible disadvan- 
tage of obtaining a singular £ matrix (to be pointed out that 
(304.l) requires,, the root of a positive number !). In the construc- 
tion selected, L and L' are always non singular ; furthermore, if 
A is the dominant diagonal, 7 and f • are equivalent. 

C C 

Finally, it may be observed tnat if the construction of ^ or 
leads to an allowable dimension in the main core of the computer, 
it is impossible to construct L for very large systems without aux- 
ilieary disks. Nevertheless, these external transfers to the main 
center are required ONI.Y ONCE during the phase of factorization. 


For practical applications, having a size of a main core which 
is not to be exceeded, it is worthwhile to choose the constant C so 






r* < 


a given percentage d/lOO of non zero elements of L q or of are 010111-/97 
orized, Therefore, since j S |qo we may define L , /in0 and L' 
as follows t ’ ' d /100 

For a given constant C, let us define k, and K' in (326) (327) 

C t. 


K c = L. * 0} 


( 326 ) 


K' - ((i,j )\l[. * 0 } 


(327) 


if n(Kp) and 
K c (resp.K^) * 


n (K^) designate respectively the number 
then the relationships between the sets 




of elements of 
^ L d/ 1 00 ,L d/j oo^ 


L d/1 00 * L C 

vithC 

so 

that 

nd^) = n(K)d/l 00 

(328) 

?, 7, 

L d/ 100 * L C 

withe 

so 

that 

°(Kq) = n(K)d/l 00. 

(329) 


By analogy to remarks 


If> d=l 00, 
If d=0, ; 


M 00/1 00 


& L' 


L 100/ 100 


are identi- l 
cal to 


correspond to the Meijerink- 
Van der Vorst type Choleski 
incomplete factorizations. 


It should be pointed out that there is another ^V 7 construct- 
ion, which is interesting theoretically, even though in 3-D appli- 
cations it leads to excessive d/lOO percentages. This construction 
is valid only for matrices using the finite elements method. 


If t? designates a standard triangulation of domain is a set 

of adjacent polyhedrals T, composed of (Mj_) N nodes. 

The complementary Ky of K may be expressed then (330) 


\ ■ {<i,j)|M.,M . cT for at T oi <£ } # (330) 

J least one 


From Ky 
construction 


it is possible 


of T rn 


to define 


in 


(33l) Kyy serving in the 


T 


*v" {(i ’ J>, 3 M K tJSt *i>V T I 

M.,Mj £T 2 

for at least one couple T^ T,, of 


( 331 ) 


Vv“ (*ij I ^ij " A ij ix £ K W } 
«>• 

l . . " 0 otherwise 
U 


(332) 


¥ith such a construction, *YV is independent from the numbering 
of • Unfortunately, the case is that L_, has few zero element: . 

within its band ( 20% in 2-D, 50% in 3-D). 

Remarks : The introduction of is also motivated by the finite 
elements method. In fact, it is easy to verify that if Q c^, then 
= 0(h) where h is the average length of the sides of T£^, where- 
as if Q c , then *£j ~ 0(1), It is also necessary to eliminate the 
small elements by a test along their width relating to the diagonal 
elements and not along thair absolute width. 


11.5 Applications of Incomplete Fac torizations to Transonic Flows and 
to the Navier-Stokes Equations 


The matrices d/100’ d/I00’“W have been introduced in the 

lifting least squares mecnods on industrial applications of large 
dimension in order to treat the algorithm ENTIRELY within the main 
cora of the computer. 

Two strategies are presented and compared with respect to infor 
mati cs (c omputation time, memory space). 

S L L ’ > 

1 d/1 00* d/1 00’" TV are used uniquely as preconditioning oper 

ators in the solution of discrete Dlrichlet problems within the algo 
rithm, thus keeping the metric H*. We have only to substitute for 
the direct descent-climb LL*, a preconditioned conjugate gradient al 
gorithm £ d/ 100 of which the convergence speed depends essentially on 
the percentage d/100, Two iterative algorithms on the pressure of 
the Stokes algorithm are presented on Flow Charts 4 and 5. 


5 t" t" » 

2 ‘'d/I 00’ d/I OQ’^T' a re used as auxiliary metrics modifying 

this time the convergence speed of the least squares algorithm. The 
direct descent-climbs LL*' are substituted by the direct descent- 
climbs In this case a minimum percentage d/lOO is required to 

keep the convergence velocities at an acceptable rate. 



2 

ITERATIVE SOLUTION ON THE PRESSURE IN L (ft) 

OF THE STOKES ALGORITHM Pl/P2 (TAYLOR-HOOD ELEMENT) 
WITH (*) SOLVED BY PRECONDITIONED CONJUGATE GRADIENT ££ t 


(S_ ) min (J(p) 
TH pc L 2 (fl) 



dx I - Au “ 
1 P 




u -zc (H^Q)) N ) (*) 


In (S TH ) P '"*• Ap - $*u is coer- 
TH p cive xn 


l 2 (Q) 


3 a >0 (Aq,q). > a||q||% V qtL (12) 


FLOW CHART k 


Preconditioned Stokes Algorithm (T-H) 


6 


3WrswV*ilH “K 


FLOW CHART 4 (cont) 


ALGORITHMS (S^y) - Initial! ation 


p°€ L 2 


- Au° - t-$p° , u°-z c <H^(n) ) N 


o ± -*-o 
l - V-u 


. O O -*’0 ■ > o 

h - g ;x - u 


ll8 n H 2 2 

p iL 

' ” (H°,h n ) 


' n+1 n . , a 
P ■ P “ P- ’i 


New direction 


n+l n £ +n 
8 ■ g - P n "*X 


IU nt! ll 

I 

n I 

8 I ■ 


h-' . g "*' ♦ ,v 


ix"* 1 • ?h n+l i x"* 1 £ (h'(!))) N 


(*) N Dirichlet problems decoupled by iteration, solved 
by preconditioned gradient A e LL 11 
N = dimension of the space 







ITERATIVE SOLUTION OF THE PRESSURE TRACE IN 
H “l/2(r) OF THE STOKES ALGORITHM P1/P2 
(GLOWINSKI-PIRONNEAU BLMENT) WITH (*) (**) 
SOLVED BY THE PRECONDITIONED CONJUGATE GRADIENT 

ALGORITHJI r L A > LL t ; A = SS* 

[ss 1 


'101 


<v 


®in | » f 

Xch‘ ,/ 2 (D/rI ^ 


AX A dr 
r 

(**) 


in < E h> * aA 
3 a > 0 


3$. 


$•? , p-A e 

f-Vp A , u-z e (H^(fl)) N | 

fcu. i A /r\\ 


\ * <t> € %(Q) 
(*) 


is coercive in h“ 5 ^ 2 (D 

^AA,A> a ot|| A IJ 2 v AcH~ 1 / 2 (r) 

- 1/2 ' 


H 


where <•,.> designates the duality 
product 


H 1/2 (f) in h“ I/2 (D 


N = dimension of the space 

* K+2 Diricblet problems in Qt solved by preconditioned 

conjugate gradient --t 

XjLi 

*'■’ t 

** doscent— c 1 inib on - with A = £S 


It 


Og 


ft 




ALGORITHM (E^ 
Initialisation 


»°'*h 


Ap° ■ , p°-X° cH^$) 

Au° - ^p°-"f , "u°-z e(H‘0B)) N 

o 

■ ^»"u° , $° cH 1 ^) 

o 

r° - r'g° 


o o fo -*o 
h « r ; X ■ u 


(h n t8 n ) 

“ " (AhV) 

A”* 1 - A" - p h n 
n 


New direction 


n+l 

g 

- g n - P Ah" 
n 


r n+I 

- X-' g"*' 


v" - 

(s n *',r" +l > 


T * 

<g".r”) 


h n *' 

- r n+l * Y h" 
n 


A *1+1 

A? 

n n+l .n+l 

• 0 , p -h c 

»>> 

*+n+ 1 

AX 

n + l ^n+1 ...1 

■ ?P . X„ e <H 0 

tf)) N 

AO"*' 

■ V*x , $ c 


Ah"*' 

3A n+l 

„ W j 

an >r 









Numerical experiences of these two strategies applied to 2-D, 

3-D transonic flo s, on the one hand, and to the Stokes algorithm 
(E. ) 2-D, 3-D, expanded in chapter 10, from the Navler-Stokes equa- 
tions, on the other hand, ore presented in chapter 12, 

12. - NUMERICAL EXPERIENCES 

12.1, Data Processing Aspects 

The numerical simulations presented below have been applied on 
the IBM 370“168 computer. 

In the case of the approximations P k , k*2, the various integrals 
involved in the derivation of nonlinear systems with finite dimension 
discrete transonic equation (t) - discrete Navier-Stokos equations 
(NS) are computed EXACTLY with FORMAC (A. LAPLACE (22)). 

For example | (T) requires the implementation of a polygonal with 
degree 3 (333/, whereas the (NS) convection terms require the inte- 
gration of polynomials with degree 5 (334) 



/L k (2L k -l) k-1,2,3 

*4L.L. k *4 ,5,6 i*j i,j-1,2,3 


(333) (T) 


(334) (NS) 


The expression of (333) (334) as a function of area coordinates 
(L^) together with their derivatives (refer to O.C, ZIENKIEWICZ (24)) 
the standard relationships (335) (336) following dimension 2 or 3 of 
the space. 



a 8 y 

L| L 2 L 3 


dr 


a ! 8 !y ! 2 
(a+y+y+2) ! 


a+8+Y * 5 , tt.fi.Y s 0 


(T) 



L i - I . 


(335) 



“AV 

'l L 2 L 3 L 4 


di 


a jB !$ »y ; 6 
(a+p+Y+S+3) 



L. 

l 


I 


(336) 


oh-B+y+S s 5 , a,8,Y * 0 

The various triangulations used --. 1 'e generated (case 2-D) auto- 
matically by the MODULEE techniques (3 r >). The large number of solu- 
tions of tho discrete Diriehlot problems justifies iho choice of a 
Choleski band or Cholcsk i-profi 1 e type direct method (35). 


l±02 


It is obvious that ths factorization phass of tho Dirichlot ma- 
trices shall always bo performed ONCE AND FOR ALL prior to the iter- /lOk 
stive process* The matrices are aotyfted entirely in the main core in 
the case of simple 2-D tests, whereas in most applications in indus- 
try 2-D/3-I>» their memorization requires ONCE AND FOR ALLdata trans- 
fers with the use of auxiliary disks. For more details , MODULEF (35) 
may be consulted. 

Finally | mention should be made of the preliminary phase of re- 
numbering the triangulation nodes t"hfor reducing the band widths of 
the Dirichlet matrices, by the CUTHILL-MCKEE algori hms (36). 

12.2. Calculations of Transonic Flows 

12.2.0. Characteristics of a Transonic Calculation 


12.2. 0.1. The_Outguts 

/ \ fff 

For each case of calculation (difference of potential) * 

incidence), we have acces, in the form of plottings, to the flow 

analysis 

-either in the fluid by the Machs distribution (337) on element* 
or the iso-Machs 


M 2 - 


2 

(y-h)I 


ifrl 2 


; $ 


( 337 ) 


-or on the bodies by the suface distribution of pressures Kp 
(intrados-extrados in the case of an airfoil profile) 



/'-£f N 2 \ Y/Y '' 


038 ) 


Remarks t 

l) The pressure and the Mach depend on the gradient of the po- 
tential. In the case of the approximation PI, the velocity is 
constant on each triangle. The Mach and the pressure on the profile 
are from two ADJACENT triangles. In the case of the approximation 
P2, the speed ( j is linear. We may therefore represent the Mach 
and the pressure on the profile by a linear variation on the bar of 
the DJACENT triangles, but a discontinuity at the inter-bars may be 
observed. ( Figure 30 ). 





•discontinuity of the pressure 
and of the Mach depending on 
(Ti,T 2 , T 3 ) 

A 

xcontinuity of the pressure and 
of the Mach depending only on 


2) The location of the shock (numerical) depends on the approx ' 
imation. 

In PI. a shock is located necessaritly at the inter-elements 
(Figure 3 1 )» whereas in P 2 it may be taken into account inside an 
element (Figure 32 ) 



12.2.0.2. Finite Elements Pl/Finite Elements P2 Comparisons 


For a same domain and a same case 

example for a flow around a circle), we 
tri angulations of Figure 33 ( PI ) and 3^ 
the schemes. 



h 


of computation ( = .45 for 

have tested the effect of the 
( P2 ) on the convergence of th 

1 



h' = h/2 


Figure 33 




Figure 34 - Triangulation P 2 


Bringing to mind the terminology "PI iso P2" : it is an approx- 
imation composed of the same degrees of freedon as the triangulat- 
tion P2, each triangle P2 gives k sub- triangles PI by joining the 
middles of the sides. 

The convergence of the schemes of optimal control formulations 
with rerrulation . penalty or artificial viscosity is verified during 
N iterations of control in the form of plottings on which are shown 

-the evolution of the cost function (C®, C*. . ...C**) . 

-the evolution of the gradient (G® , G* , . . .G N ] ; G" a (g^g* 1 )*'^ 

-the determination of the circulation (joukowski condition) 

-the determination of the physical shock (supersonic-subsonic 
domain) 

-the local . action of the penalty terms to prevent the devel- 
opment of shock decompression. 

12.2.0.3. F'JNITE^ELEKENTS^FXNITE^DIFFERENCES^Comgarisons 

The unconservative and conservative codes of A. JAMESON have 
served as reference for numerica. tests on the NACA 0012 airfoil 
and the KORN airfoil. 

It has proven to be instructive to compare locally the shock 
INTENSITY and LOCATION in lifting and non lifting cases between the 
two conservative codes (Finite Elements + Penalty) and (Finite Ele- 
ments+ Artificial Viscosity) of the optimal control and of the two 
JAMESON codes (Consorvatife Finite Differences) and ( conservative 
Finite Differences) at 1 50 degrees of freedom (on the airfoil pro- 
file) and iso case of computation ( and identical incidence). 

Moreover j, the difficulty of treating the Joukowski condition 
in finite elements (p + =p“* measured in AVERAGE at trailing edge in PI f 
exactly on airfoil profile in ^2) was able to be disconnnected from 
comparisons (Finite Elements PI - Finite Differences ) by calculating 
with iso CZ (CZ : aerodynamic reaction of airfoil). Most of the re- 
sults which follow have already been presented either in GP^lB (37) » 






or in contrac c LABORIA/IRIA/DRET (38). 

12.2.1. The Jonverging-Di verging Pipe h°i 

The potential $ is given at the pipe inlet and outlet , whereas 
the condition of tangency _ y is applied on the sides. 

3n 

The domain of the flow is quantified in 384 TRIANGLES on figure 
35 for an approximation PI - rough card-index or P2 (resp 1536 in the 
case PI ISO P2). 

The number of corresponding nodes was 221 (resp 825) for a lin- 
ear approximation (quadratic resp. or PI ISO P2°. 

Figures 36 and 37 give a comparison without condition of entropy 
and with condition of entropy treated by REGULATION with P = . l(respp 
= .2 , u. = .1) of local Machs on axis ( 3 ) and the side (4) of the 
pipe. “ H 

4o iterations (resp 60 ) were required to obtain the convergence 
of the conjugate gradient algorithm thereby requiring 1.30 ran of pro- 
cess (resp ?mn). 

Figures 38 and 39 show a plotting of the iso-machs resulting 
from P2 measurement in the regions (subsonic -supersonic) and (super- 
sonic - subsonic) of the flow with shocks. 

The agreement of the two approximations may be verified. 

12.2.2. The circle /l 1 3 

The NON LIFTING flow around a disk has a double numerical value: 
the equal distribution of the points of quantification on the circle 
due to a constant curve and of the compression and decompression 
shocks with equal intensity located symmetricallly . We have select- 
ed a case of transonic calculation Mo = .45. 

For this problem the boundary conditions are the NEUMANN type 

= U,„ .n at infinityt = 0 on obstacle). 

3n 3n 

Numerical considerations require the substitution of a bound do- 
main for the infinite domain with r, sufficiently far from the ob- 
stacle in the following sense : if j is the chord of the obstacle, 
the distance of l from the obstacle is equalt to about 4 or ^ times 


The domain is divided into 3456 TRIANGLES (resp 834) corres- 
ponding to 1813 NODES for one linear approximation (resp. quadratic). 

The condition of entropy was treated by PENALTY and the conver- 
gence of the algorithm is obtained in 50 iterations (resp, 60 ) cor- 
responding to 4 mn of process (resp. 8 inn). 
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Figure 40 shows a comparison PI ISO P2/P2 of pressures (Kp) 
calculated on the ADJANCENT triangles to the circle* 

Figure 4l shows in an iso-mach form the case in P2 of the shock 
on a single element ADJACENT to the circle* 

It may be observed that there Is a strong shock intensity for 
the computation case y obtained by the two codes* 

12.2.3. The NACA 0012 Airfoil Section (Profile) 

A rough triangulation brought about by a WINSLOW algorithm (39) 
(reap* fine) (60 points on the airfoil) with enlargement near the 
obstacle , is given on figure 43. It is composed of 1080 triangles /l 1 7 

(reap. 4380) and 600 nodes (reap* 2280 ) . “ 

12,2.3.1. The symmetrical non lifting cu^e (without JOUKOVSKI condi- 
tion 


Two test cases have been calculated t 

■ (M.n “ *3 ; INC » 0*) "non scitf" case 

(2) - (>U - .85; INC » 0*) • "stiff", case 


Figures 44 through 48 relate to (l) 

On figures 44, 46, 47 we have plotted the distribution of pres- 
sures on the airfoil profile. 


The results of figure 44 (resp. 45) correspond to a treatment of 
the condition of entropy with PENALTY (resp ARTIFICIAL Viscosity + 
REGULATION) ( y . .00i5 ; K - .4) (resp v - .05 ; y - .00000*) 


In the two cases, the convergence of the conjugate gradient al- 
gorithm was obtained in 40 iterations corresponding to 3*5 mn of pro- 
cess. One may notice the clearness of the shock obtained with Pen- 
alty. 


Figure 46 compares the solution obtained’ by PENALTY in PI on 
the fine triangulation with the ones derived from the conservative 
and non conservative codes of JAMESON in finite differences . 

A comparison in the sense of approximation PI ISO P2/P2 is made 
on figure 47 with the PENALTY ( PI ISO P 2 :y*. 5 ;K*. 0 ; 

P2 : Uj * . 1 ; K - .4 ; y, * .01). 

One may take note of the shock case in P2 on a single element 
ADJACENT to the airfoil profile together with the rocompression after 
the shock, which marks the conservative form of the equations. The 
iso Machs near the airfoil profile derived from computations PI and 
P2 with the entropy-penalty condition have been plotted on figure 48 
and give an idea of the location of the shock and of its intensity 
in the fluid. 
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The effect of the artificial viscosity is represented on 
figures 58 “ 59 vith tho local Machs noar tho trailing odgo ( 58 ) 
and in tho shock rogion (59)* 

A Fintos Elomonts comparison is givon ( on a rough and fins 
tri angulation ) on figure 60 . "" 

It may bo observed that the quality of tho compression shock is 
restored by the PENALTY at the supersonic - subsonic passage. 

Finally, a result (P2 - PENALTY) with predictor PI ISO P2 of 
artificial viscosity type gives a good result on figure <>1 • The 
supersonic zone of the two calculations PI and P2 in the fluid , 
in the vicinity of the profile defining the position and the intensity 
of the shock is represented on figure 62 ; it may be observed that the 
shock is taken into account in P2 on a single element adjacent to the 
profile. 

The interpolation problem Pl/P2 makes it possible to give to 
code P2 a good predictor PI and is presented in (4o). 
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Figures 49 through 52 relate to ( 2 ), 

The PENALTY has been used on figure 49 with U*' • The conver- 
gence is obtained after 60 iterations corresponding to a process 
time of 4 mn. 


The local effect of these terms of PENALTY during the itera- 
tions is shown at the bottom of the decompression shock on figure 50 . 

It may be pointed out that at the end of the computation! the 
constraints remain active and this brings to light the unstable na- 
ture of the solution. 


A comparison in the sense of the approximation PI ISO P2/P2 
(u * l./u* l and ^2 * .01) plotted on figure 51 . The location and 
intensity of the shock on the airfoil are shown by the iso-Machs of 
computations PI and P2 on figure 52 . 

12.2.3*2. - The Lifting Case ( With JOUKOVSKI Condition) 


Two test cases have been calculated : 
<3> (M„ * .6 ; INC = 6 °) 


- Small supersonic zone f but strong inten- 
sity decompression shock very near the com** 
pression shock. 


- Large Supersonic Zone. 

(4 ) (Mo, * .78; INC » I®) 

Figures 53 through 56 relate to (3)* 

Figure 53 compares the pressions on the airfoil with the JAME- 
SON finite differences non conservative and conservative method with 
the pressures obtained in PI with ARTIFICIAL Viscosity + REGULATION 
(v* .005, h * .00001) on a rough triangulation. The local Machs 
in the shock region at the extrados of the airfoil are shown on fig- 
ure 54. A comparison of the supersonic zones in the form of iso 
machs P1/P2 shows a good agreement between the two Eipproximations on 
figure 55. 

A PI Finite Elements comparison ( PENALTY- VISCOSITY (ARTIFICIAL) 
on figure 5 6 brings to light the good behavior of the code with PEN- 
ALTY which at the same time in a very narrow zone, restores the phy- 
sical shock and resists the high intensity decompression shock. 

Figures 57 through 6l relate to (4). 

Figure 57 compares the JAMESON finite differences conservative 
and non conservative solution with the solution obtained in PI with 
ARTIFICIAL VISCOSITY + REGULATION (v « .005, \i •= 5.10~ 6 ) on a fine tri- 
angulation. 20 iterations on the JOUKOWSKI condition have been per- 
formed, representing 80 optimal control iterations for a process time 
of 15 mn . 
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12. 2. *4. The KORN Airfoil Section ^Profile ) 



The Korn airfoil section is a non symmetrical section designed 
to produce a transonic flow without shock if .75 and INC = 0°, 

Since the flow is not symmetrical, the JOUKOVS condition is applied 
at the trailing edge. 

The domain of calculation surrounding the section has been divi- 
ded into 2880 triangles (resp 1362) for a piece-wise linear approxi- 
mation (resp. quadratic) with 1560 NODES of which 120 on the section. 

The triangulation with detail near the section is provided on 
figure 63 . 

Figure Gh shows an effect of the triangulation (rough and fine 
on the location and intensity of the shock for the test case M,, a .75 
INC = 0° with artificial viscosity + Regulation v = 005 ; 

y * .00005). 

It may be observed that the shock intensity decreases with h, 
quantification step. 

Comparisons (Finite differences, JAMESON conservative scheme) - 
( PI finite elements (rough triahgulation) )- (P2 finite elements) - 
are presented on figure 65 . The condition of entropy was treated by 
PENALTY. 60 iterations for a process time of 30 mn are required to 
obtain the convergence in case P2. 

Another case of computation with iso CZ (M - 75 ; INC = O'. 1) 

disconnecting in order to treat the J0UK0VSKI condition demonstrates 
the agreement of finite elements + artificial viscosity with conser- 
vative JAMESON finite differences on figure 66 . 

A second test case has been performed m = .75 and INC = .5 /l4l 

OH "" wmmmmmm 

A comparison Finite differences - Finite elements with PENALTY 
( u = . 1) at 150 degrees of freedom is presented on figure 67. Att- 
ention shall be brought to the compression shock clearness of the 
solution with penalty . 

Finally, the conservative case H.,= .75 and INC = .5, calcula- 

ted either by the JAMESON Finite differences, or by the PI Finite 
elements with artificial viscosity ( v = .008, V* - . 00005 )duct of very 
similar solutions on figure 68. 
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12,2,5. The Multibody (NOZZLE ♦ AIRFOIL SECTION 1 


/148 

The lifting transonic flow around an industrial configuration of 
multlbodlea has bssn calculated and compared for approximations PI, 

PI iso P2 and P2, 

The triangulation around the MUL 2 for a linear approximation 
(resp, quadratic) consists of 2936 elements (resp, 734) corresponding 
to 1 553 Nodes* The matrix factorized by Cholevski is composed of 
200 610 coefficients (reap, 256 276) whereas the number of non zero 
coefficients of the DIRICHLET matrix is 10533 (resp 17725)* Details 
of the rough triangulation of the nozzle and of the slot is given on 
figure 69 * The Joukovaki condition is applied to the trailing edges 
of the nozzle and of the airfoil section* 

The condition of entropv is treated by REGULATION* 

The test case H m «.5l TNC»10° calculated on MUL 2 ? FINITE /149 

ELEMENTS PI ♦ REGULATION < - " •*> (resp. P2 u,- .5 & V 2 » .05 
required 80 control iterations corresponding to a process time of 15 
mn (resp, 23 mn). 

Figures 70, 71 , 72 show the surface Mach distribution on the 
nozzle (l) and the airfoil section (2) for rough triangulations Pl y 
fine PI iso P2, and P2* One may see the presence of a shock at the 
extrados of the airfoil section (2). 

Details near the mulibodies of the local Machs in the fluid in 
the form of the Mach number (Pi) or iso-Mach (P2) shows the g od oper 
operation of the nozzle, the passage at M«1 at the neck on figures 73 
74, 75 (downstream from the slot) and the satisfactory Joukowski con- 
dition on the nozzle (l) at the subsonic limit ! * .95). 

The determination of the circulations during the iterations de- 
pending on the approximation selected is shown on figure 76, whereas 
the evolution of the cost function and of the gradient of the criter- 
ion depending on the approximation selected are compared on figures 
77 and ?G. 

It may be pointed out on figure 78 the "periodic” discontinuity 
of the gradient corresponding to the calculation of a now circulation 
(joukowski condition) and requires a restoration of the conjugate 
gradient algorithm in the sense of POWELL (4l), 

12.2.6, The BI-NACA Multibodv AIRFOIL SECTION ♦ AIRFOIL SECTION 

The interest of a transonic calculation around a (BT-NAC) con- 
figuration lies in the mixed naturo of the simnl taneeusly interna.l- 
extornal flow. In fact, the internal domain ('u) made up by the ex- 
trados of the lower airfoil section ( 2 ) and intrudos of the upper 
airfoil section (l) is the converging-diverging pipo type, vhoreas 
the one (flj) formed by the intrudos of (2) nnu b, the extrados (l) 

(Figure 79*1) represents an external flow around a body. 
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The poeaiblity of several shocks appearing simultaneously and having /i6Q 
different Intensities is therefore a fundamental numerical test for 
industrial applications 2-D 3-D composed of several shocks* 

The triangulation around the BINAC is made by the MODULEF tech- 
nique and consists of 3298 elements corresponding to 1739 nodes* The 
number of non zero coefficients of the Dirichlet matrix is 11.800, 
the one factorized by Choleski has 147*117 coefficients. Details of 
the triangulation near the 2 airfoil sections showing the internal 
domain is given on figure 79*2* The Joukowski condition is applied 
to the trailing edges BF1 and BF2 of the 2 airfoil sections 1-2* 

Two cases of computation taking up 2500K of di ole precision 
memory I) (M^ * .6, INC - 0°) (non lifting) 

2) (M u . - .6, INC - 6°) (lifting) 

in finite elements PI with PENALTY are presented and have required 80 
iterations corresponding to a process time of 20 mn* 

Figures 80-81 show the surface distribution of the pressures on 
airfoil section 1 and airfoil section 2* It may be observed that 
there is a perfect symmetry of results o) the pressure intrados of (l 
(l) is mixed with the pressure extrados of (2) and vice versa* Case 
1 has only one shock inside the domain pipe type, whereas in 

case 2) a second shock is placed extrados of (2), in the external 
domain . 

Details near the two airfoil sections of the local Machs on 
figures 82-83 in the fluid, in the iso-Mach form show a good opera- 
tion of the internal domain and the satisfaction of the Joukow- 
ski condition* The penalty prevents simultaneously the formation of 
two decompression shocks* 


12*2*7* The Converging-Diverging 3-D Pipe 


Z±£7 


This is an ajustment test case of code 3-D* The appearance of 
compression shocks is verified in the diverging part of the pipe, 
as the formation of decompression shock was prohibited by the penalty 
of the condition of entropy. 

As in came 2-D, a difference of potential is applied at the in- 
lit and outlet of the pipe which is sufficiently high to obtain a 
case of transonic operation* On the sides, the tangency conditioner* 
* 0) of homogenous Neumann standard type are implicity applied, 1 " 
The domain of the flow is quantified into 1920 tetrahedrons on fig- 
ure 134. and is composed of 24 sections* 40 iterations lead to con- 
vergence of the algorithm in 3 mn of process time. 


on 


On figures 85 through 90 may be seon the evolution of the Mach 
numbers, constant on each tetrahedron, on several fronts ad j cent to 
the sections located in the converging zones (without shock) and di- 
verging zones (with shock) of the pipe. One may note the satisfact- 
tion of the entropy condition in a region near the pipe axis. 
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12. 2. 8 . The 3-D Air Inlet 


This is an industrial application. Figure 84 shows a detailing 
of the TETRAHEDRONI Z ATI ON (see J.G. NAVES (52)) near the air inlet 
in a vertioal plane, whereas on figure 91 there is shown the geometry 
of the air inlet together with one part of the tetrahedronlzation 
(fronts of the tetrahedrons attached to the air inlet) used for a 
piece** wise linear approximation of the potential. 

The external and internal domains of the air Inlet are made up 
of 25664 tetrahedrons corresponding to 5732 Nodes . To give an idea 
of the complexity of the problem, one may observe that the Cholevski 
matrix l(A-LL c ) (of the discrete Dirichlet operator) is made up of 
about 2 million coefficients and that its factorization requires 15 
mn of processAime. »» 

REGULATION' has teen used to treat the condition of entropy. 

The computation test case^® • .3 ; \ 0cOr ■ -55 ; INC * 6* » 

DERAP s 0°) has required 40 iterations corresponding to 60 mn of 
process time. 

Figure 92 gives the Machs internal and external surface distri- 
bution on the tetrahedrons • DJACENT (in the direction of one front) 
at the air inlet and shows the narrow supersonic band on the upper 
external part of the air inlet. 

The long computer usage time, due to inputs-outputs of the 
factorized matrix L, is the reason for the incomplete numerical fac“ 
torization tests £ presented . in paragraph 12.4. 

12.3.0. Characteristics of an Incompressible Viscous Calculation 

Tba.Iasuia 


In the velocity-pressure formulation, a Navier-Stokes calcula- 
tion required of boundary conditions on ^ and sometimes on p. Three 
situ tions are encountered in the applications. 


1 . 

2 . 

3. 


Dirichlet conditions on the entire boundary T , u| 

3u. ^ 

* 3n'r 


Neumann conditions on - one part of 
4 Dirichelt condition on the pressure 

Mixed conditions on the comoonents 
of velocity : 


L - 0 


3u . 


but such that 
7«u « 0 in 2. 


•*i ; 3n T. 


- 0 j ^ 


u«n dF ■ 


s s 

0* constraint required by the condition of 

compressibility, 

. -* 

* • an external cal- 


As the equations are without dimensions l u r.,! 
culation around an obstacle (airfoil section - air inlet) requires 
the assumption of an incidence and of the Reynolds number p m J_ 
with u fluid viscosity, c v 


v 







The Reynolds Is reduced to a characteristic lenght L, which in 
the example shall always be the unit (cavity lxl, diameter circle 
del | chord profile {»i ,i air Inlet deviation h=l). 

As the velocity and pressure approximations may vary (figure 93) 

i l) pressure PI - velocity PI 
2) pressure PI - velocity PI SIO P2 
3) pressure PI - velocity P2 

one computation requires the simultaneous presence of two card i nd— 
ices in the computer corresponding to two triangulations CC*,€*l,/ 9 ) 
for example. The discrete Dirichlet operator A is therefore n * 1 '- 
constructed twice, depending on whether it is applied to the pressure 
or to the velocity • Furthermore, it may be observed, in 

the unsteady case, that it depends on the time step Et and on the 
Reynolds number Re, since in this case the metric of the generalized 
Stokes algorithm is expressed 

" <f» C (kId-vA)$ . (339) /1 79 

The two tri angulations h and h/2 are therefore numbered twice 
by the Cuthill-MacKee algorithm in order .to obtain band widths m^ 
and m 2 at minimum 

Ibfi-QutBUta 

The (velocity-pressure) formulation permits direct acces to the 
fields of vel oc ities (l) and pressures (2) and to the vorticitv in- 
tensity ( 3 ) 7 a u constant on each element if ^ is PI, piece-wise lin- 
ear when y is P2, The streamlines (4) are obtained by solving a Dir- 
ichlet problem (340) at a riven field of velocity 


- Aij. ■ V a u (ft) 

<H r " 8 (D 


(340) 


The visualizations of magnitudes (l) (2) (3) (4) in the form of plot- 
tings at various time cycles i\t make it possible to follow the evo- 
lution oi the flow in limo (origine of eddies, appearance of spoara- 
ted zone, alternating emission of eddies in the fluid, corresponding 
pressure fluctuation on the bodios). 

The plotting of the iso-streaml ineo , tho iso-pressures and the 
iso-vorticities is ensured by tho T 7 L\C0 modulus ( refer to MAUItOCCO - 

IN TKll LID (42)). 






As the values ^M J N and ^MAX are determined after solving ( 340 ), 

N desired values of 0 ? i«! ,N, with possible cubical concentra- 
tions on pax-licuiar .(V“ 0 > values ••• (change of sign mark- 

ing the eddies in the iluid) are marked geometrically (x,y) on the 
bars of triangulations rs with a linear connection from one 

bar to aiother on each h ^h/2 * element. 


The convergence of the schemes of approximation is verified 
during N control iterations by plotting 

-the evolution of criterion (C° t C* , , , . ,C^) 

-the evolution of gradient (G° ,G* , , . . ,G^) ; G^ = (g^»€f^)^'^ 

-the values of constraint 
-► 

The control w is initializd following the applications either 
by the solution of the Stokes algorithm, or by the idealized fluid 
solution. In external flows, the Stokes solution proves to be a 
poor predictor. 

In the unsteady case, the sequence of optimal control problems f 

is initialized at the solution of the preceding time cycle, each ! 

problem requiring a few control iterations if the time steps are not ; 

too large. 


Finally, industrial applications require numerically high lam- 
inar Re ynol ds (R e = 1000), a climb in Reynolds by a parabolic law of 

f\ lk”l 0 type shown on figure 94 makes it possible to 

simulate in a wind tunnel the transitory 
phase of determining the solution by Rey- 


1 00 


(k=l ,10) 


nolds calculations. 


For each value v i of the viscosity, the matrix is not re- 

constructed (which would be a penalty in computer time), but is sub- 
stituted by an equivalent modification of the velocity boundary con- 
ditions expressed in figure 94. 


Most of the following results are shown in R. GLOWINSKI-B. MAN- 
TEL - J. PERIAUX-O. PIRONNEAU (43), in IRIA/LABORIA-DRET (l9), AMD/ 
BA-DRET (44). 


12,3.1. The 2-D Test Cavity 


The Stokes flow in a cavity (lx l) was tested to verify the 
error estimates of schemes o ( h a 1 of BERCOVIER- PIRONNEAU (45/ for 
three approximations (Pl/Pl) (Pl/Pl ISO P2) (P1/P2) of the (pressure- 
veloc7.y) formulation. 


The characteristics of the 3 triangulations 


studied 




corresponding respectively to the values hg 
are defined in (341 ) 
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M^T^*3PSS»af*?£- m - 


^ Q 

K nodes* 256 elements} 


C M 


j, “ {221 nodes* 400 elements} 


(340 


^ P 

. = {841 nodes, 1600 elements} 


The calculation of Qt def ining the convergence of the scheme 


is obtained to satisfy the constraint -*■ _ n evaluated numerically 
in (342) and (3*3). V * U ” ° 


DIVGL0 


L f |W|** 

T c*? h >a 


DIV MAX - sup (|?.IL) 

Ttf. T 

n 


( 342 ) /1 81 

(343) 


For each approximation, a is defined for the possible couples 


(G,M), (o,p), (M,P) by the formulas (344) ( 3^3) (346) 


a(G,M) = 


a(G t P) = 


a(M,P) 


Log 

DIVGLO 

(G) 

DIVGLO 

(M) 

Log 1 .25 


Log 

DIVGLO 

(G) 

DIVGLO 

(P) 

Log 2.5 


Log 

e — - - 

DIVGLO 

(M) 

DIVGLC 

(P) 


(344) 

(345) 


Log 2. 


(3M) 


For data on the edge of the cavity (u = 16 x^(l-x)^, v=o), 
we have plotted on figures 95 and 96 a Log-Log scale, the slope of 
a of the straight line Log Divglo = cl Log h. characterizing the 
scheme 0(h a ) depending on the approximation chosen for the GLOWINSKI- 
PIRONNEAU Stokes algorithm and the optimal control Na' r-Stokes 
method at Re = 100, after 30 control iterations, Scht O(h) is 
verified approximately for the case Pl/Pl ISO P2 and O(h-) for the 
case P1/P2. For more details, (46) may be consulted. 
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12.3.2. 


The 2— D Conduit With Sudden Enlargment 


The characteristics of triangulations 
about by MODULEF ( 35 ) » are given on figure 


(f & 

h 

97. 



brought 


/1 86 


It may bo observed that there is a concentration of elements in 
the recirculation zone. The calculation domain (6x » Sy) . the 
boundary conditions (u = zi(y) } v=0) and the Reynolds number are 
proposed by A.G, HUTTON (47), Two cases Re = 100, Re = 191 are ob- 
tained by making the unsteady code steady Pl/Pl ISO P2 in 180 iter- 
ations corresponding to one time step At and requiring 3b. of 

process time. 


Superposing the 
ure 98) shows a good 


Si, = jx 


enlarge- 

ment 


streamlines with those of I;heHUTTON code (fig- 
agreement along the length of the blister (if 

and h designates the height 
of the enlargement 


x point of 
connection 


6JL = 6 x h a Re * ICO, 6i. ■ 8xh Re » 191). 


The appearance and the developement of the separated zone throug 
through various time cycles At a t Reynolds 100 are shown on figure 


On figures 99 through 102, the field of pressures and stream- 
lines of the flow at the two Reynolds numbers under consideration 
may be compared. 


12.3.3. The Alternating Eddies Behind the Circle /1 93 

The The triangulation T»^(resp. ^h/2^ is composed of 144 elements and 
84 nodes (resp. 576 triangles and 312 nodes), the solution ) 

looked for is composed f 708 degrees of freedom. “ ^ 


At Reynolds 50, the Navier-Stokes solution is steady, as the | 

streamlines show on figure 103, after 40 time cycles. Nevertheless, | 

with this Reynolds, the Stokes solution is already a poor predictor j 

on figure 104 at time cycle 1, 

At Reynolds numbers above 80, the steady solutions of Navier- 
Stokes equations being unstable, we consider the unsteady case, as * 

the flow is initialized at t=0 by the incompressible idealized flow. | 

Since the approximation keeps the symmetry and the triangulations ! 

tf h/2 are also symmetr’ ;al, the solution (u h 10 *PhV s shown on fig- i 

ure 105 (a) corresponding to K=10 (t= 10 At) is symmetrical and must 

therefore be perturbed at a point of fluid not round on the axis. Ac- j 

cordingly, we may observe behind the circle the formation of a Karman 1 

path. The results presented on figure 105 (a)-(f) correspond to ) 

Reynolds Re = 200 and are obtained by an implicit Orank-Nicholson j 

type scheme with a time step At ~ .1. The process computation time j 

is about 1 hour. We have verified the good agreement of the results ' 

obtained by FORTIN-THMASSET (48) by using a different unconform mix- 
ed finite elements method. 
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Fib res 106-107-108-109-110 show a flow Re ■ 200, with Nemann /19T 
condition behind the circle, which is inadequately perturbed to pro- 
voke alternating eddies in 50 tine cycles, 

12. 3, *4. Separated Flow At Extrados of Airfoil Section In Incidence /202 

We are taking into consideration an unsteady flow around an 
airfoil with Reynolds 200, placed at 30° incidence. 

The calculation domain is substituted by a triangulated bound 
domain by M0DULEF(^. : 412 triangles, 221 nodes j £ h/2 l 1648 tri- 
angles, 854 nodes). The solution looked for (t? >p ) is composed of 
1929 degrees of freedom. The quantification in time is accomplished 
by a completely implicit Gear scheme with two steps, with one time 
step At - .1. The predictor *0 is the solution of the incompress- 
ble idealized fluid. u h 

80 time cycles (corresponding to a period of 8 seconds) have 
required 90 mn of process time and u core space of 1500 k octets. 

The number of control iterations per cycle of time is 4. 

The velocity distribution and streamlines on figures 111 (a) -(f) 
and 112 (a;-(f) show the formation of eddie ixtrados of the airfoil 
which alternately expand and escape in the xxuid to be finally absor- 
bed by the downstream boundary conditions. 

12.3.5.1. The Air Inlet In Incidence /207 

The mixed flow around inside an idealized air inlet with high 
incidence is a typical example of a separated viscous flow. Iu a 
first phase, the air inlet is placed at an incidence of 30 ° as is 
shown on figure 113 . 
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There are two types of boundary conditions verified by the /207 
velocity t 

-Dlrichlet on *« ^ea * u “ (q) 

-mixed Dirichlet-Neumann on 


u*l 

3v 

3ST* 


Thus we define the velocity satisfying the constraint 

_ « where n designates the external per- 
il* n * 0 pendicular to r. 

‘“ r EA ur - wr s 

The domain p. is triangulated by the MODULEF techniques ( 35 ) • Th 

The triangulations £ and ^ the characteristics of which are given 

. « n/2 

on figures 114-115 , are relatively rough , but on the other hand, 
they cannot sustain a large Reynolds number (Res! GO, Re reduced to 
h, distance of the 2 airfoil sections 1-2). 


12.3.5.2. Solution of the Stokes Algorithm 


/208 


In a first phase, we have compared from the point of view of 
informatics (calculation time) and of theory (accuracy of the scheme) 
the solution of the Stokes algorithm either by mixed formulation (u,0) 
FLOWINSKI— ^IRONNEAU, or by the TAYLOR-HOOD formulation (u,0> • The 
first approach relates the the numerical solution of (E ) expanded 
in 10.6.5.3. by a conjugate gradient iterative method on the pressure 
trace \ on r , whereas in the second one, the conjugate gradient al- 
gorithm is used on pressure p in ft , described in R. GLOWTNSKI-O, PI- 
RONNEAU ( 49 ) . 

The two algorithms converge for a same approximation P1/P2 
toward a pressure distribution in ft which is very similar, on 
figures 116-121 after satisfaction fot the stop test on g” i 
g n : (6 n .gV /2 < e, in 30 iterations. ( c . , 1D _ 6 ) 

The conditioning S occurring in the solution (3^2) (3^3) is 

2 n 

taken in L , optimal choice in the TAYLOR-HOOD approach, since 


V"* 1 ■ V" - p V “ ith Vh ■ "s' 

V"*’ ‘ V" * P V" With V ■ '"“hq 


(3 42 ) 

(343) 


j 


-I /2 

but not in the GLOWINSKT-PIRONXEAU one, sinco ' ' 1! i • ) • Wo 
thorofore expect to improve the convergence speed of ( 3 ^ 2 ). 
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12.3.5.3. Comparisons of Codes Pl/Pl ISO P2 and P1/P2 


/217 


Comparisons of the two codes are based on the following case i 

-unsteady Navier-Stokes flows with Dirichlet or Neumann condi- 
tion downstream. 

Two cases have been calculated if i = 30° , Re * 100. The time 
step selected is At - the number of time cycles selected is 100 f 

the number of control iterations at each At is 6 . The process compu- 
tation time is about 100' in the P1/P2 case, 55' in the Pl/Pl ISO P2 
case. 


It may be stated that on the whole the numerical simulation of 
the flow obtained by one or the other code is very similar. 

Figures (122) (123) (124) show through the means of streamlines 
at Re = 100, the appearance, the development and the discharge of 
large structures on the upper external part of the air inlet and in 
the internal part, the formation of a quasi-steady eddy, which re- 
mains attached to the lower side. 

Since the domain of calculation is voluntarily selected to be 
small, the boundary conditions downstream interfere considerably 
with the entire flow as soon as the ejected eddies reach the down- 
stream boundary, which is shown by the gobal flow at time cycle 100 
(velocities, streamlines and pressure of figures 125-127 (resp. 128 - 
130) For Dirichlet type conditions (resp. of Neumann type). 

Interpretation of the results confirms the choice of Neumann 
type downstream boundary conditions for larger Reynolds. 

It is interesting to observe the .numerical operation of the two 
codes by following the evolution of values of criteria and gradients 
through time cycles and within one of them. It may be observed that 
when the Reynolds number increases, it takes longer for the conver- 
gence of the optimal control problem to be obtained (3 to 4 itera- 
tions for Re * 50, whereas 6 to 8 iterations on the average for Re = 
100 ). 


Figures 1 31 through 1 33 show the evolution of the criterion and /21 8 
of the gradient within a time cycle without much alteration in the 
flow. The following 134 through 1 36 figures relate to a time cycle 
( 75 ) close the the emission of a new eddy. 

It may be observed that code P1/P2 absorbs •'better" the altera- 
tion in configuration, whereas code Pl/Pl ISO P? shows more resis- 
tance (jump of criteria and of gradients) and requires more itera- 
tions to control the new fluid state. 


20G 


p l / P| lio P t 




" Time Cylcle w Time Cycle 

‘ Reynolds ico. Reynolds 

STREAMLINES 

NEUMAN DOWNSTREAM 


1 ± 




Figure 122 




/ P* l«e P % 

to Time cycle 
'^•Reynolds 

STREAMLINES 



y o Time cycle 
l0C) - Reynolds 



7 5 Time cycle 
,30> Reynolds 


33 Time cycle 
ioo Reynolds 

DIRICHLET DOWNSTREAM 



so Time cycle 

ico. > 

* Reynolds 



ioo Time cycle 
l00, Reynolds 


Figure 123 



















TIME CYCLE 





CRITERIA OR 
GRADIENT 



Pf / P 2 /228 - 

CRITERIA AND GRADIENT CURVES 

AT THE BEGINNING OF TIME ITERATION 



Figure 131 

Re s 100 



CRITERION 


TIME CYCLE 










CRITERION 
OR GRADIENT 




: 


P < / P 1 ISO P* 


CRITERIA AND GRDTENT-eDR'"- - 
AT THE BEGINNING OF TIME ITE! 


TERM 


Gradient *. 

CRITERIA* 


Gnid/ent I 


CRITERION 


no-.IOG 

\ 

4 


Figure 123. I 


/ 


X-:/ 


TIME CYCLE 





CRITERIA AND GRADIENT CURVES 
AT TIME CYCLE. 7 5 . . 






ioo\ 


FIG. 1 34 

Gradient 


CRITERION 



I 


GRADIENT 

& CRITERION 


20 Qj 


n 

\ 

\ 

\ 

r 

\ 

\ 

\ 

\ 

\ 

• : \ 

\ 

\ 


Pi /Pi ISO Pi 

i ♦ * • i , » I 

CRIJTERljfli ANP'‘GrCADI^f,-r CURVES' 
AT. SEIMd-CldLE -7:5 i L .. 



t : 

l. L. : 

I ! ; 

1 

l ; ; 

I 

t ” i 

l 

l • • 

I 

• i • : 

\ 

• \ 

\ : 

\ 

. V ... 

\ 

\ 

- \ 


fig. i: 


J ' Grad/enf 


CRITERION 


HC 1 GO 




i 



TIME CYCLE 




Plnally t It may be observed on figure 1 37 that, in code 
P1P1 ISO P2 t the Neumann type condition downstream facilitates 
the measurement of the alteration in configuration (smaller jumps 
of criteria and of gradients). 

The comparison of the process computation times between the 
two codes y brings to light a ratio of 2 in favor of Pl/Pl ISO P2, 
this figure is directly related to the amount of calculation for 
the creation of various second members according to the approxima- 
tion P , k*»l or 2 through time cycles and especially to the amount 
of quantified Laplacien Choleski coefficients (as the band width m2 
of P2 is about 2 times higher than for band width m2 of Pl/Pl ISO P2. 
In the case under consideration m2 s 129, ml ■ 68, the corresponding 
core space is 987 K for the case P2 and 540 for the case PI. 

Ve shall see that, given the Reynolds range considered in in- 
dustrial applications, the compromise Pl/Pl ISO P2 is a sensible 
choice. 


12,3.5.4. The Industrial Configuration i«40*. Re » 250 /% 

The operation of the air inlet, proposed by ONERA (refer to H. 
VERLE (503) around/in which is simulated the separated flow, has 
been studied experimentally in the form of visualizations with Rey- 
nolds 3 ( |q 4 The case computed (Re ■ 250) j i*40 # ) is composed of 
6893 degrees of freedom. TriangulationsCT, , , -created automatically 

by MODULEF (35) are shown on figures 1 38-1 39. ^Phe density of the 
nodes near the air inlet is shown on enlargments 140-141, The large 
amounts of factorized discrete Dirichlet matrices requires the use of 
auxiliary disks with a Choleski "shyline" FLIP-FLOP method escribed 
in MODULEF (35). 

Due to the high incidence, a parabolical flow c « . g inside the 
air inlet (percentage of j u^)- applied in order to prevent a poss- 
ible blocking and to suck vno eddies formed at the air suction inlet. 

100 time cycles calculated with a time step At « .05 have required 
several hours of process time. 

Figures 142 (a)-(f) (velocities), 143 (a)-(f) (iso-pressures) 
show the formation, the development and the ejection of several ed- 
dies inside and outside the air inlet. It may be seen on figure 144 
(f), which represents the streamlines, the existance of 5 eddies with 
alternating signs, cf which 2 are inside the air inlet spreading a- 
long the entire height and sliding slowly toward the aspirator ! It 
may also be observed that the streamlines in front of the air inlet 
are drawing closer together, which will effect the quality of the 
approximation, (density of nodes) the more the Reynolds is higher. 

One may have a better idea of the complexity of the flow by 
looking on figure 1 46 at the superposing of the time cylce 100 of 
142-f, 1 43-f , 1 44-f , 
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12,3.6. The 3-D Sphere 


l^L 

The incompressible viscous fluid flow around a sphere with dia- 
meter 1 proves to be an interesting: informatics test to check the 
satisfactory operation of a 3-D Navier-Stokes code from the symmetry 
properties of the flow y the obstacle and the tetrahedron formation. 

The conditions applied to the boundaries are the Dirichlet type 





0 

0 

0 


Informatics problems due to the 3-D and to the analysis of re- 
sults on this example are immediatly sufficient. The domain of com- 
putation 8 is formed into a tetrahedron containing,^ 2 4 elements and 
154 nodes in PI (G h ), » decomposing: in Pl/ISO P2 0^/2' into 4992 

elements as on figure 147 and into 970 nodes (reaching thus 2000 de- 
grees of freedom the solution (u^p^) obtained by the optimal control. 

Minimization of the band width proves to be an essential pre- 
liminary step if we want to work with factorized Dirichlet matrices 
having a size acceptable in the main core, requiring 60 * process and 
1900 K of core space. 


Since the time step is At ■ .1, 40 time cycles at Re » 100 is 
sufficient to induce behind the sphere a separated zone shown on fig- 
ure 148, 


Visualization of the return velocities is shown from the side 
and globally by hachuring the tetrahedrons, of which the component 
of the velocity ^ is negative. 

12,3,7, Swept-back Ving at Large Incidence /248 


In this industrial example, we have taken into consideration the 
3-D flow of an incompressible viscous fluid at Re = 200, around a 
complete left-right idealized wing, placed at 30 ° incidence. 


The triangulation % consists of 2060 tetrahedrons and 560 nodes 
Due to the importance of the factorized Dirichlet A matrix (A ■ LL C , 

A constructed from an approxin n P^), 74562 coefficients, we are 
focusing in a first phase on a j.,.near approximation of the velocity 
on constructed like a). We are assuming that there are enough 

node- m L for us to solve (E. ), 

f l h 


The calculation (70* process) consists of 40 time cycles, with 
the time step being At s ,| , the number of control iterations at 












Fig. 148 






each cycle being The use of auxiliary disks for solutions AX ■ b MS 
brings us to consider two different computer times t one time t. of 
process for the computation volume itself and one time tg machine 
space due to external trans'ers to the main core t„ ** nt with 

I -ns ID, highly dependent on the informatics environment at the 
moment of computations. 

The solutions (*J,p) at various time cycles are registered on disk 
to be analyzed after the computation. Visualir.. - : ;ions make it poss- 
ible to identify the separated zones which are obtained in the fol- 
lowing manner. 

1. Several angles are plotted (views from the front, side, rear 
from above, below, in perspective) at various time cycles, the set of 
tetrahedrons T c in which the component u of velocity V » (u,v,w) 
is negative . The support of the entire wing is represented by a 
plotting with a different color, making it possible to locate the se- 
parated zones and to evaluate the intensity of them (figure 149 (a) 

(b) (c)). 

2, From a separated zone, we can plot the lines upon which the 
vorticity is applied (vorticity tube lines) to visualize the eddy 
intensity (A, MARROCCO ( 51 ) ) . 

Depending on the starting point (end of the wing, for example), /gfro 
we may represent, in the separated zone, the complex path of the 
fluid particles. 


Various views of the three dimensional eddies are shown on fig- 
ures 150 (a)-(d), I 5 I (a) (c) corresponding to two integrations with 
different initial conditions. 


On figure 1 50 (a)-(d), we are focusing on eddies which escape 
at the end 'f the wing (left or right), whereas on figure 151 (a)- 
(c), we are placed initially in a less turbulent separated zone. 

It may be stated that the two separated zones interact, since the 
integration of the vorticities from the wing-right provides traject- 
ories leading to the separated zone of the wing-left via the socket. 


The numerical integration of the lines is obtained by the fol- 

‘■•T of XZ 


lowing process : given a point Z„ of the separated zone 0 0 

we calculate the vorticity cj_ *= * a u_ constant in T 0 , the velocity u 


l! 


being P 


1 ' 


Since we,., are looking i 

* n 


the geometrical intersection Z, 


-* with fronts/;- '> 

W u i 1 i*l .4 

F,_ found gives a new tetrahedron T 


of T and a point Z 0 at Z 


1 * 


The front' 


. . (close to T # in the dir- 
ection of F ), Ke calculate the new vorticity ■+ of element T. and so 
forth... k w, 1 
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12.4 Data Processing Efficiency for Functional Least Squares Algo - 
rithms Using Auxiliary Operators or Metrics 

12.4*0. Characteristics of a Preconditioned Computation 

The examples of 3-D transonic and incompressible viscous flows 
have demonstrated the need for auxiliary disks which serve to memor- 
ize the factorized Dirichlet matrices L (= 5. 10° coefficients). These 
matrices are read numerous times during the uescent-climbs of a so- 
lution LL*X a B and result in excessive memory use time (5 times more 
than process time). 

It may be recalled that one optimal control iteration requires 
5 Dirichlet solutions in the transonic case and 35 Dirichlet solu- 
tions in the Navier-Stokes case. 

The objective of these examples is to show that preconditioning 
operators L{j/ioo» constructed in chapter 11, make it possible to 
solve entirely in main core a problem which initially exceeds the 
computer capacity. Vue present two ways to use conditioning operators 
in an optimal control problem. 

l) The matrix is kept in the penalty (344) or B plays the 
role of the discrete Laplacien 


min 

4>eR r 




BE 


BE = R($)} 

(*) 


(344) 


B ~ L L* 

but, d/100 d/100 d/100 is used as auxiliary operator of Laplacien 

in the sense of 0. Axelson to solve (*). In this case, the conver- 
gence speed of the algorithm is not slowed down, 

2) The metric H^- is approached in formulation (345) by g g 
plays, then, the role of the auxiliary metric . 5 


min N , (E^E 
$eR-' 


BE = R($)} 
(**) 


( 345 ) 


but in this case, it shall be fair to choose percentages of B d/100 

such that d/I 00 > d°/1 00& so the algorithm does not slow down excess- 
ively. It nay be pointed out, on the other hand, that (*#) has„an 
extremely fast solution : a descent-climb of one operator ^d/ 1 00 Lq 0 
representing, for example, 20 of the Laplacien if d=20, The choice 
of d in case 2 is the better compromise between (345) and the pos- 
sible size of the computer. Examples of ^ f are shown on figure 
152. Attention shall be brought to the d/l0 J proximity of 
non zero coefficients of j- d-100 kept to those of A. 
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12.4.1 a Auxiliary Operators and Metrics in Transonic 


J25k 

12*4,1,1, 2-D Laplaciens Preconditioned by {'£ c 

In a first phase, it is worthwhile to test a conjugate gradient 
algorithm to solve the problem (346) by the finite elements method. 


miiv. F$) ( 34 6 ) 

$eR 


in which t \ = LL C is introduced as an auxiliary operator of the La- 
placien operator A in the sense of 0 AXELSSON (32), 

is constructed from the factorization L (11267 coefficients 
of A and L^/ioO represents various percentages of ~ constructed in 
accordance with the procedure described in 11. 

We have plotted on figure 1 53 the number of iterations required 
to solve (346) with a specified accuracy c =.10-6 , by using L ( i/|C0 

and l' d'/lM constructed in (324) (325) for different d's. We may 
note the interest of the interval ( 5 %, 25?0 for memory decrease, and 
compare the convergence velocity with other auxiliary operators such 
as the Van der Vorst operator L^, , which does not require factori- 

zation L or still constructed oy keeping only the coefficients 
very close to L and representing a small percentage (205c) in 2-D. At 
both ends -of the curve, we find the solution of (346) in one itera- 
tion for ^ 1 00 / 1 00 and the standard conjugate gradient, without pre- 
conditioning. 

On figure 1 54 . we have superposed two curves d/100 d/100 

as a function of the number of iterations with and constructed 
in (324), but from two different renumberings of the Cuthill-McKee 
algorithm j L contains 11267 non zero coefficients and L' 13569* The 
agreement of the two curves may be verified when working with iso- 
percentages on the two auxiliary operators. 

12.4.1.2. 3-D Laplacien Preconditioned by 

The solution of (345) has been also found on an industrial con- 
figuration with 5328 degrees of freedom, of which the Choleski ma- 
trix L contains 1.5 D° coefficients and could not be held in the main 
store. 


Figure 155 ~describes the number of reasonable iterations 
when operators ^/lOO’ d/100 < 20/100 aro used in the main stor» 

of the computer. We jnay note the number of excessive iteratioi 
\ 1 4 6 2 ) of the standard conjuage gradient when (346) must be sol 
several hundreds of times. 
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12.4.1.3. 


Transonic Optimal Control 2-D With Metric 
Operator ™ 


and Auxiliary 


The approach 12.4.1.1 is used to solve the state equation (*) of 

(344) 


BE ■= R«>) 


(*) 


by preconditioning the conjugate gradient algorithm by d/I 00' We 
shall point out the safety of the algorithm J344) which converges in 
N iterations regardless, of the conditioning ^/ioo selected to solve 
(*). We have shown on figure 1 56 the process computation timo to 
perform- a transonic computation on a NACA 0012 at " *8 ; i*0°)by 
using L d/10Qf°r several values of d. We shall bring our attention 
to the interest of the pointss of the curve in the interval 5%$ 2596 
producing about the same process times as those using high percen- 
tages. The optimal control formulation using the standard conju- 
gate gradient as Lsplacien algorithm is very costly. The curve sta- 
bility is kept by working on another numbering of the triangulation 

12.4.1.4. Transonic Optimal Control 2-D With Auxiliary Metric LL fc . 


When the metric attached to the solution of the transonic oper- 
ator is perturbed in the sense of (3^5) » N iterations required 

for a transonic computation may increase if the metric b = Ll^is too 
weak (percentages too low of d/100 comparing the initial metric H^- 
with the metric L^), 


Figure 1 5 7 shows for various choices of d/lOO the evolution of 
the error g, taken in the good standard g t BC» committed to solve the 
equation R(c£>) = 0 in the functional space H“^-, as a function of the 
control iterations. 

When the metric Hi/lOO is acceptable in the sense of the con- 
vergence, the solutions of (3^5) prove to be faster and more econ- 
ical in store than the standard solution. 


It may be observed that the Van der Vorst 
operator used as auxiliary metric to solve an optimal control^ pro- 
blem via(345) is inadequate. On the other hand, the metric L^ v « 
composed of very close coefficients and representing in 2-D about 
20 <jc of the coefficients of L, appears to be an acceptable auxiliary 
metric on figure 157 . 


2.H2 


l^L 


i 


The quality of the transonic solution as a function of the 
various auxiliary metrics (various percentages d/lOO, Van der Vorst 
after 80 control iterations is represented by shock restoration on 
the airfoil section ( on figure 158 . It may be concluded that 15# 
is the minimum allowable percentage for an auxiliary metric. It 
still represents a considerable gain in memory for industrial 
applications. 

12.4.1.5. 3-P Transonic Optimal Control With Metric H^-and Auxiliary /256 
Operator ~~t . 

The solution of (344) using the preconditioning ^d/100 of (#) 
has been tested on an industrial type air inlet configuration com- 
posed of 1.5 10“ Choleski coefficients and 5328 degrees of freedom, 

at M « .8. 

00 

The curve of figure 1 59 represents the .process time of N»10 
control iterations for percentages d/100 ofL^yjQQ entirely in the 
main store. We may note the vertical slope of the curve as soon as 
d/100>5%, f expressed by, the constant number of iterations required 
to solve (*)- AS LONG AS L d/100 i s in rHE MAIN CORE. The point ob- 
tained with L ioo/I00 and an auxiliary disk depend on the working con- 
figuration of the computer at the moment the computation is performed. 
Fluctuating usage times may be obtained for the same calculation at 
various phases. 


12.4.1 .6. 3-D Transonic Optimal Control With Auxiliary Metric £Jt . j 

The same industrial configuration has been tested by solving 
(.344) via (345). The error evolution for various auxiliary metrics 
B d/100 is shown on figure 1 60 during the control iterations. It 
may be seen that is the minimum metric leading to an allow- 

able error curve compared to reference Bjqq^jqo* 

Since the Van der Vorst metric ®vdV is too far from the fact- 
orized L of the Laplacien, it is poorly suited for the solution of 
(345) and leads to an insufficient convergence velocity. 

It may be concluded, after examining figure 160, that d/lOO ■= 

20 £ ^ is an auxiliary metric making it possible to treat (345) entirely 
in the main core and to ensure the convergence of the preconditioned 
algorithm with a sufficient safety margin. 
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12.4.2. Navler-Stokes C .so Around /in an Air Inlet (2-D) a :-d Around 
A 1-D Wing 

12.4.2.1. The Stokes Algorithms (T-H) and (G-P) 

12.4.2.1.1. Preconditioning of the Lan^acien In the Iterative 
Stokes Algor ithm ( T-H ) 

We havo introduced in the iterative algorithm of Flow Chart 4 
(see Chapter ll), a preconditioning ‘|L t ln the solutions in 2-D and 
3-D of the Dlreichlet problems. Figures 161 and 1 62 show the evolu- 
tion of calcluation time for solving the Stokes algorithm (T-H) with 
accuracy e - . 1 0—6 - given on the pressure , for various precondition- 
ing percentages L d/ 100 * It may be pointed out that the optimal work- 
ing zone, hachurated on the figures, the economy S/ 100 < d/ 1 00 < 20/1 00 
of memory (* 90~) does not penalize at all the computer process time ! 

The 2-D example (resp. 3-D on the shpere) was initially composed of 
9342 (reap. 149734 ) Choleskl coefficients on the air inlet for the 
factorized matrix L. Algorithm 4 does not call for preconditioning 
on the pressure, since the conditioning L 2 ln the Taylor-Hood ap 
proach is optimal . 

12.4.2 .1.2. Preconditioning of the Laplacien and ss c of the 

Pressure Trace i t the Iterative Stokes Algorithm (G-P) 

We have introduced in the iterative algorithm of Flow Chart 5 
(see Chapter 11 ), first, a preconditioning a > LL 1 ^ in^the solution 
of Dirich 1 at problems, second, a preconditioning a » SS C on the pres- 
sure trace , since the conditioning L 2 in the Glowinski-Pironneau ap- 
proach is not optimal. 

Figures 163. 164 show in 2-D and 3-D the eveolution of calcula- 
tion time for solving the Stokes algorithm (G-P) with accuracy c 10-6 
given on the pressure trace, for various preconditioning percentages 

Lj/lOO “ nd S d/I00' • Sln<! * " atril A i ‘ < =°"' plet ' *4/100 iS ° b ‘ 
tained by a test, absolute in 2-D, and relative in 3-D, on the amount 
of coefficients of factorized A. 

Figure J 63 shows the optimal working zone, hachurated, corres- 
ponding to L 24 /,oo< d< 6/100 and S 34/|00 -* Tt ma y be observed that the 

preconditioning \2 inadequate. Figure 164 shows the fast 

decline in computation" time in 3-D, as soon as percentage of g 
which is too small, is used. 

If a comparison is made of the calculation time of the two ap- 
proaches (T-H) and (G-P), it comes to light that it is better to 
work on the pressure trace (factor 3 to 4). 

In any case, the numerical tests shown on figures 160 through /268 
1b3 clearly show' that the procondi tioning problem of a Dirichlet op- 
erator in H *s perfectly solved., whereas the problem of a t ra c e 
operator on I* remains open. 
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12. 4,2. 2. Optimal Control (,?-D)(3-P) Navier Stokes metric H^ - - /2 73 

Auxiliary Operators ££t 


In industrial applications (2-D) (air-inlet) and 3-D (wing), 
informatics memory problms are due to the storage of the Dirichlet 
operator ( a id-vA) , on the one hand and of the trace operator^ -*■ aA «* 

on the other hand. 


the 



r 


/T\ ~y. 

An alternative u ' w in order to gain is memory space is proposed 
for solving the Navier-Stokes equations via Flow Chart 1. 

CD Apply the direct algorithm of the Stokes algorithm (G-P)(Flow 
Chart 3) with preconditioning l\} to solve the sequence of Dirichlet 
problems. In this case, we must construct upstream of the optimal 
control loop the trace operato r A aA , symmetrical but complete , with 
the use of auxiliary operator » Chen factorize it (A = SS 1 ) and 

store S (flow chart 2). The importance of S may require auxiliary 
memories for direct solutions of Ej, ; aA = b with the auxiliary oper- 
ator 7 being completely stored in tne main memory. 

L d/100 

(J Apply the iterative algorithm of the Sotkes algorithm (G-P) 

(Flow Chart 5) with proconditioning s °l ve the sequence of Dir- 

ichlet problems. In this case, a preconditioning sS t the trace A 
operator is necessary in order for the time required for solving, 
compared with the direct method, is still competitive. Nevertheless, 
making the choice remains delicatel Two auxiliary trace operators g 
are suggested. 

^ 2 

2.1. We use S^.(A) * A N.*dr, conditioning L , restricted to the 

boundary node supports of figure 2U, With this choice, operator A 
is never constructed . It may be observed that - is sparse , its 
memorization presents no problem. 

.2.2. We use a percentage ^d/100 of the complete matrix A = SS* after 
constructing the latter upstream. For various percentages d/100 re- 
lating to the relative -.value of coefficients (S±a) (j > i) , w& obtain 
conditioning operators ^d/lOO of which the efficiency is measured 
a posteriori by the convergence velocity. 

The use of auxiliary operators LL and SS C in a Navier-Stokes / 27 1 * 
algorithm is presented on figures 1 65 (2-D) and 166 (3-D). We have 
set in a* ssa the process time for treating tho Navior-Stokes com- 
pletely 4... the main memory in 10 iterations with a small Reynolds 
number (Re - 50) os. a function of percentages d/lOO and d'/lCX) of 
operators and S. Attention may be brought to the fast increase 

in calculation time for percentages d'/lOO of ^ which are insuffi- 
cient ( d ' i 50) . On the other hand, for a given S, the interest of 
operators “d/100 (5sd<20) may be pointed out, as they have very 
little effect on tho process time, while representing a gain in mem- 
ory space of about 9 0 %! 
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CONCLUSION 


The quality of the numerical results of this study confirm that / 277 
the functional least squares methods coupled with preconditioned 
conjugate gradient algorithms proves to be a tool which is particu- 
larly suited for multiple industrial configurations. The possibility 
of treating correctly the boundary conditions of any complex geometry 
by a finite elements method, gives to the codes obtained from the me- 
thod presented, a flexibility which is indispensable to the three di- 
mensional aerodynamics of today and of the future (optimum design) . 

In the case of transonic flows, it appears that Lagrange P ap- 
proximation by conform finite elements (resp. mixed) of the related 
optimal control problem, including the condition of entropy treated 
by penalty (resp. artificial viscosity), is of sufficient accuracy, 
after comparison with results derived from the A, JAMESON finite dif- 
ferences codes. 

With respect ot the incompressible viscous fluid flows, the com- 
plexity of the Navier-Stokes equations suggests the use of quantifi- 
cation schemes of lower order, for velocity and P^ for pressure. 

The convergence, however, is ensured only if the tr : 5 angulation of the 
domain used for the velocity is twice as fine as the one required for 
the pressure. 

In the two flow' families considered, the approximation by mixed 
finite elements (artificial viscosity in idealized fluid, Stokes al- 

f orithm in viscous fluid), as presented in P.G. CIARLET-P.A. RAVIART 
54), R. GLOWINSKI (55), GLOWINSKI-LIONS-TREMOLIERES ( 56 ) and J.M. 

THOMAS (57) for the biharmonic problem and more recently in FORTIN- 
THOMASSET (48) by the Navier-Stokes equations, remains a very impor- 
tant point. 

Sophisticated codes, obtained from the optimal control-Stokes 
algorithm combination, and the convergence of which is ensured by the 
absolutely stable CRANK-NICHOLSON .implicit schemes, while being per- 
haps more costly in machine time and memory usage, are easier to use 
in industry (no convergence parameters to set ! ) than the tradition- 
al codes requiring domains of reduced stability. 

The numerical simulation of three dimensional separated large 
structures, the dimension and location of which play a fundamental 
role in aerodynamics with large incidence (interaction of eddies em- 
itted by several bodies, life-time of eddies in the air inlets) is a 
demonstration of feasibility of the optimal control tool, which is 
indispensable in the subsequent phase of combining Navier Stokes with 
turbulence models . In any case, calculations with a large Reynolds 
number is still prohibitive, if not impossible, with the size of com- 
puters currently available (sequential organization of computations), 
the memory capacity of which proves quickly to be inadequate for as- 
sociated quantification ( 1 0 () calculation points for 3-U applications 
is not an excessive number ! ). 
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The incomplete factorization methods presented in the nonlinear /278 

context (solution of the Dirlchlet problem several hundreds of times! 
brings a gain in memory space of the order of a factor 10. Introdu- 
ced in the conjugate gradient algorithms coupled with optimal control 
in the form of auxiliary operators (preconditioning - ll c the Dir- 
ichlet problem LL C $ - F) or auxiliary metrics (minimization in jp 1 of 
F(i) ■ C) , they make it possible to solve entirely in the main memory 
3-D configurations taken from the two flow families , and this is ac- 
complished in acceptable machine times. 

They represent, however, only an intermediary stage, if compared 
with the possibilities of parallel calculators of tomorrow. 
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